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ABSTRACT 


With  the  aim  of  developing  a  fast  and  accurate 
computer  code  for  predicting  the  aerodynamic  forces 
needed  for  a  flutter  analysis,  we  review  some  basic 
concepts  in  computational  transonics.  The  unsteady 
transonic  flow  past  airfoils  in  rigid  body  motion  is 
adequately  described  by  the  potential  flow  equation  as 
long  as  the  boundary  layer  remains  attached.   The  two 
dimensional  unsteady  transonic  potential  flow  equation 
in  quasilinear  form  with  first  order  radiation  boundary 
conditions  is  solved  by  an  alternating  direction  implicit 
scheme  in  an  airfoil  attached  sheared  parabolic  coordinate 
system.   Numerical  experiments  show  that  the  scheme  is 
very  stable  and  is  able  to  resolve  the  highly  nonlinear 
transonic  effects  for  flutter  analysis  within  the  context 
of  an  inviscid  theory. 
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I.   INTRODUCTION 

We  begin  with  a  survey  of   the  behavior  of  flows 
past  conventional  airfoils.  Then,  we  introduce  the 
unsteady  transonic  problem  in  flutter  analysis.  Finally, 
we  discuss  the  mathematical  difficulties  of  solving 
such  a  problem. 


1 .    Flow  Past  Airfoils 

We  begin  our  discussion  with  a  brief  survey  of  the 
behavior  of  flows  past  airfoils;  when  a  conventional 
synmetric  airfoil  accelerates  from  subsonic  speed   to 
supersonic  speed   the  flow  pattern  usually  develops  in 
the  manner  shown  in  Figure  11.   As  the  flight  speed 
of  the  airfoil  reaches  the  critical  speed,   the  local 
flow  speed  equals  the  local   sound  speed.   Beyond  the 
critical  speed,  a  supersonic  region   appears  on  the 
airfoil  which  is  usually  terminated  by  a  nearly  normal 
shock  through  which  the  flow  speed  jumps  from  super- 
sonic  to  subsonic.   With  a   further  increase  in  the 
flight  speed,   the  shock  moves   aft   and  the  size  of 
the  supersonic  region   and  the  shock   strength   both 
increase.   If   the  pressure   jump  through   the  shock 
is  sufficiently  large,   separation  of   the  boundary 
layer  occurs.   This  shock  induced  separation  starts  when 
the  local  Mach  number,  the  ratio  of  local  flow  and  sound 
speeds,  just  upstream  of  the  shock  is  about  1.2  5  to  1.3. 
When  the  boundary  layer  downstream  of  the  shock  separates, 
the   nature  of  the  flow   around   the   airfoil   changes 


completely  and  often  turbulent  flow  phenomena,  such  as 
buffet  or  buzz,   start  to  occur. 

The  other  important  flight  parameter  is  the  angle  of 
attack,  the  angle  between  the  flight  direction  and  the 
airfoil  chord.   The  effect  of  changing  the  angle  of  attack 
of  a  conventional  symmetric  airfoil  at  a  given  super- 
critical speed  is  shown  in  Figure  9  .   When  the  angle  of 
attack  is  increased,  the  speed  over  the  upper  surface 
increases,  and  the  shock  strength  and  the  supersonic  region 
on  the  upper  surface  both  increase. 

The  flow  patterns  for  a  modern  supercritical  airfoil 
acceleration  in  speed  or  angle  are  usually  similar  to  the 
patterns  shown  in  Figures  12  and  10,  respectively.   The 
supersonic  zone  in  these  cases  may  consist  of  several 
pieces. 
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2.    Engineering  Considerations 

An  aircraft  under  certain  circumstances  may  experience 
vibrations  of  an  unstable  nature.   This  phenomenon,  called 
flutter   in  aeroelasticity ,   is  governed  by  the  interaction 
of  the  elastic  and  inertial  forces  of  structure  with  the 
aerodynamic   forces  generated  by  the  motion  of  the  vehicle. 
These  forces  interact  in  such  a  way  that  the  vibrating 
structure  extracts  energy  from  the  passing  flow.   This  may 
lead  to  a  progressive  increase  in  the  amplitude  of  vibra- 
tion and  may  cause  structural  damage  and  loss  of  control  of 
the  vehicle. 

For  a  given  vehicle,  the  aerodynamic  forces  increase 
rapidly  with  the  flight  speed  while  the  elastic  and 
inertial  forces  remain  essentially  unchanged.  There  is  a 
critical  flight  speed  called  the  flutter  speed,  above  which 
flutter  occurs.   The  requirement  that  a  flight  vehicle  be 
free  of  flutter   over  the  entire  flight  range,   which  may 
include  subsonic,  transonic,  supersonic  and  hypersonic  speeds, 
is  one  of  the  most  crucial  factors  in  the  design  and  construc- 
tion of  flight  vehicles.   The  vibration  characteristics  of 
the  vehicle  at  zero  speed   can  be  determined  quite  accurately 
by  numerical  methods  or  ground  vibration  tests  [44] .   Thus 
flutter  analysis  depends  mainly  on  the  knowledge  of  the  aero- 
dynamic forces.   In  subsonic  and  supersonic  flight,  aerodynamic 
forces  can  be  predicted  reasonably  well   by   current 
methods  based  on  linear  theory.   For  transonic  flight. 


nonlinear  effects  make  the  evaluation  of  the  transient 
aerodynamic  forces  considerably  more  difficult.   This  has 
concerned  the  flutter  analyst  since  the  beginning  of 
transonic  flight.  The  transonic  regime  with  its  mixed 
subsonic-supersonic  flow  patterns,  usually  containing 
shock  waves,  is  the  most  critical  regime  for  the  deter- 
mination of  the  flutter  boundary.   A  typical  flutter 
boundary  with  transonic  dip  is  depicted  in  Figure   1. 
The  flight  speed  may  exceed  the  flutter  speed  in  the 
transonic  region. 


speed 


1.4 
Mach  number 


Fig.   1.   Typical  Flutter  Speed  vs.  Mach  Number  Curves 
of  a  Flight  Vehicle. 


Currently,  supercritical  v/ings  make  it  possible  to  cruise 
at  transonic  speeds  with  low  drag.   This  leads  to  a  renewed 
interest   in   transonic   flutter   analysis.   In  this  paper 
we  consider  inviscid  unsteady  transonic  potential  flow 
past  airfoils  in  rigid  body  motion  with  the  aim  of  provid- 
ing  a  method  of  predicting  the  aerodynamic  forces  needed 
for  a  flutter  analysis. 


3 .    Mathematical  Problem 

In  mathematical  terms  we  find  solutions  to  a  partial 
differential  equation  that  describes  flow  outside  a  wing 
section  which  is  in  rigid  body  motion.  There  are  several 
difficulties  in  this  problem: 

1.  The  equation  is  nonlinear, 

2.  The  physical  time  direction  is  not  the  time-like 
direction  of  the  equation  when  the  flow  is  supersonic, 

3.  Shock  waves  occur,  and 

4.  The  body  surface  is  moving  in  time,  which  is  equivalent 
to  saying  that  there  is  an  essential  singularity  at 
infinity  in  the  airfoil  attached  reference  frame. 
While  much  progress  has  been  made  in  the  mathematical 

theory  of  transonic  flow,  many  basic  questions  remain  open. 
For  example,  even  for  the  small  disturbance  equation,  one 
of  the  simplest  nonlinear  mathematical  models,  it  has  not 
been  shown  that  the  problem  is  well  posed  in  a  suitable  class 
of  weak  solutions.   The  linear  theory  is  deficient   in 
predicting  important  features  of  transonic  flow  outside 
airfoils  in  low  reduced  frequency  motion  [29] . 

At  present,  a  very  effective  way  to  study  unsteady 
transonic  flow  is  to  obtain  approximate  solutions  by  compu- 
tational methods.  We  overcome  the  first  difficulty  by  the 
use  of  finite  difference  methods.   This  allows  the  solution 
to  be  advanced  in  time  by  solving  a  sequence  of  linear 
equations  which  approximate  the  nonlinear  equation  if  the 


time  step  is  small.   The  second  difficulty,  as  well  as  the 
third,  is   solved  by  a  type  dependent  differencing  strategy 
which  employs  central  differencing  for  all  terms  at  sub- 
sonic points   and  upwind  differencing  for  the  streamwise 
derivatives  and  central  differencing  for  the  transversal 
derivatives  at  supersonic  points.   Shocks  are  captured 
automatically.   The  fourth  difficulty  is  solved  by  using  a 
coordinate  system  in  which  the  airfoil  is  fixed.   The 
far  field  then  has  an  essential  singularity  that  can  in  turn 
be  treated  by  introducing  radiation  boundary  conditions  at 
the  artificial  boundaries   which  are  a  finite  distance  away 
from  the  body. 


4.    Plan  of  Work 

The  plan  of  this  work  is  as  follows;   In  Section  II, 
several  flow  models  derived  from  the  conservation  laws  of 
fluid  dynamics  and  the  proper  constitutional  hypothesis  are 
reviewed  in  decreasing  order  of  complexity.   We  begin  with 
the  Navier-Stokes  equations  and  step  down  to  Euler  equations, 
potential  flow  equations,  small  disturbance  equation,   and 
low  frequency   small  disturbance  equation  .   We  discuss  the 
proper  boundary  conditions  and  related  concepts  in  each 
flow  model.   We  also  review  some  basic  numerical  concepts 
and  discuss  a  splitting  technique  for  constructing  stable 
implicit  schemes. 

In  Section  III,  we  restrict  our  attention  to  the  poten- 
tial flow  equation  in  quasilinear  form.   We  study  the 
characteristic  surfaces  of  the  equation  and  derive  the 
proper  radiation  boundary  conditions  for  the  artificial 
boundaries  of  computational  domain.    We  also  discuss 
coordinate  transformations  which  render  the  airfoil  surface 
lying  along  a  portion  of  coordinate  surface   in  the  compu- 
tational domain. 

In  Section  IV,  we  construct  a  highly  stable  alternating 
direction   scheme  for  the  potential  flow  equation  in  the  compu- 
tational  domain.   The  finite  differencing   strategy  and 
approximate  factorization  technique   are  analyzed  through 
linear  models,  convection  equation   and  wave  equation  . 
It  is  shown  that  the  scheme  is  unconditionally  stable  for 
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these  two  cases. 

In  Section  V,  we  check  the  scheme  by  calculating 
steady  flow  past  some  airfoils.  The  computational  results 
show  that  the  scheme  is  very  stable  and  there  is  no  problem 
in  calculating  sonic  flight.   Then,  we  demonstrate  our 
ability  to  calculate  unsteady  transonic  flow  past  realistic 
airfoils  in   rigid  body  motion. 

In  Section  VI,   we  present  the  conclusion  of  this  work. 

In  Appendix  A,  we  describe  a  5-diagonal  matrix  solver 
employed  in  our  scheme. 

In  Appendix  B,  we  explain  the  operation  of  the  computer 
program  and  the  glossary  of  input  parameters.   We  also 
present  the  listing  of  the  computer  program  UFL05. 
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II.   BASIC  CONCEPTS 

With  the  object  of  developing  a  fast  and  accurate 
computer  code  for  unsteady  transonic  flow  past  airfoils 
we  review  in  this  section  some  basic  mathematical  models, 
including  governing  equations  and  boundary  conditions 
for  unsteady  transonic  flows   and  some  relevant  numeri- 
cal concepts   including  the  resolution  of  the  finite 
difference  mesh  system,   the  ideas  underlying  the  split- 
ting technique   and  the  shock  capturing  technique  used 
to  construct  an  alternating  direction  implicit  scheme 
and  the  advantages  and  disadvantages  of  conservative 
and  nonconservative  difference  schemes. 
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1.    Customary  flow  Models 

In  many  aeronautical  applications   turbulent  flow  is 
observed.  The  phenomenon  of  turbulence  is  not  well  under- 
stood  and  currently  much  attention  focusses  on  finding 
useful  models  to  describe  turbulent  flow  theoretically. 
Continuous  flow  models  have  been  found  adequate  to 
describe  a  large  class  of  flows  of  practical  importance  [32] 

1.1   Navier-Stokes  Equations 

With  the  proper  constitutive  approximations,  the 
conservation  laws   of  mass,  momentum  and  energy  lead  to 
the  Navier-Stokes  equations  in  Cartesian  x,y  coordinates 
in  the  conservation  form  [32,38] 


(1) 


U^  +  F   +  G   =0 
t    X    y 


where 


U  = 


P 

pu 

pv 


F  = 


pu 


pu  +  a 


PUV  +  T 


xy 


9e 


(e+a^)u  +  Ty^v  -  K   g^ 


G  = 


pv 

puV  +  T 

2 

pv  +  a 


yx 


(^-^^y)^'-' V^'  ^If 


with 
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•3u 


9v^ 


xy 


and 


9u  ,  3vi 


3u 


X         ^9x    9y^       9x 


•9u 


9v^ 


"^'P-  ^(1^*   IJ)  -  2« 


9v 


y   '  ^9x  '   gy*    "^  9y 

in  terms  of  density  p,  pressure  P,  velocity  components 
u  and  V,   viscosity  coefficients  X    and  y ,   total  energy 
per  unit  mass  e,  specific  internal  energy  e  and  coeffi- 
cient of  heat  conductivity  k.   To  close  the  system  we 
adjoin  the  equation  of  state  p  =  p(£,p).   The  simplest 
equation  of  state  is  the  polytropic  relation  (y-law) 


P  =  (Y-l)ep  ,    Y  =  constant  , 

where  y    is  the  ratio  of  specific  heats,  equal  to   1.4 
for  air. 

The  above  system  can  be  rewritten  in  the  nondimen- 
sional  form  [42] 


(2) 
where 


U  = 


U^  +  F   +  G   =  R~-'"(R  +  S  ) 
t    X    y    e    X    y 
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pu 
pv 
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F  = 


f  0 
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G  = 
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with 


and 


T    =  (A  +  2y)u   +  Xv 
XX  X     y 

T    =  y  (u  +  V  ) 
xy      y    X 

T    =  (X  +  2vi)v   +  Au 
yy  y     X 

-1        -1     2 

R,   =  UT    +  VT  ,  +  kP   (y  -  1)    3  a 
4       XX      xy      r  X 

S;,   =  UT  „  +  VT    +  icp"-'-{y  -  l)""""  9  a^ 
4      xy     yy     r   '         y 

P   =  (Y-1)  [e  -  0.5  p(u^  +  v^)] 


where  the  local  sound  speed   a   is  given  by 
a^   =  y(Y-  1)  [e  -  0.5(u^  +  v^)  ]  , 


A  is  taken  as  -(2/3)y,  the  Stokes  hypothesis.  Note  that 
the  nondimensional  reference  quantities  are  arbitrary,  the 
Reynolds  number  R  and  the  Prandtl  number  P  used  in  equa- 
tion (2)  are  defined  in  terms  of  these  reference  values. 

Usually,  two  types  of  boundary  conditions  must  be 
specified  to  determine  flow  past  airfoils  in  motion. 

a.  The  body  surface  condition  requires  the  flow  velocity 
relative  to  the  body  be  zero  (no  slip  condition) ,  and 

b.  Appropriate  far-field  boundary  conditions  must  be 
specified  at  the  necessarily  finite  limits  of  the 
computational  domain. 
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1 . 2   Euler  Equations 

If  viscosity  and  heat  conduction  are  neglected,  the 
flow  equations (2)   are  reduced  to 


(3) 


U.  +  F   +  G   =0 
t     X     y 


and  the  equation  of  state  for  the  y-law  gas. 


(3-) 


P  =  (Y  -  1)  ep 


In  the  inviscid  flow  field,  if  there  are  surfaces 
of  discontinuity,  the  solution  of  the  differential  form 
(3)  has  to  be  interpreted  as  a  weak  solution  of  the  flow 
equation  with  proper  entropy  condition.   u(x,y, t)  is  a 
weak  solution  of  differential  equation  (3)  if 


T 


W^'U+W  -F+W  -G 
t     X     y 


dx  dy  dt  = 


W*U  dx  dy 


for  any  smooth  test  function  W(x,y,t)  which  vanishes 
for  II  (x,y)  II  large.   An  equivalent  statement  of  weak  solu- 
tions of  the  differential  form  (3)  is  that 

a.  The  differential  form  (3)  holds  in  the  smooth  region,  and 

b.  Across  any  surface  S  of  discontinuities,  the  following 
jump  condition  holds; 


n^[U]  +  n  [F]  +  n  [G]  =  0   on   S 

,.  ,-.       ;  ,.   t  X  y  .  .    .':■■.■■' 

Here  n  =  (n  ,n  ,n^)   is  a  unit  normal  vector  to   the  surface  S 
X   y   t 
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of  discontinuity  pointing  from  the  region  (1)  to  the 
region  (2) .   More  specifically,  if  q  is  the  velocity 
vector  of  the  flow  and  s  is  the  velocity  of  the  surface 
of  discontinuity,  then  the  jump  relations  derived  from 
the  conservation  laws  of  mass,  momentum  and  energy  are 


(4)  m  =  (n«q^-s)Pj^  =  (n'q2-s)P2 


(5)  mCq^-q-j^)  =  n(p^-p2) 

and 

e^    e  ^  ^ 


Equation  (5)   implies  the  following  two  equations: 

(7)  m(n'q2-n«q^)  =  Pj^  "  P2 
and 

(8)  m(nxq2  -  nxq^)  =   0 

We  can  distinguish  two  cases  with  either  m  7^  0  or 
m  =  0   across  the  surface  of  discontinuity.  In  the  first 
case,  the  tangential  velocity  component  nxq  is  continuous 
across  the  surface  which  represents  a  shock  wave;  in  the 
second  case,  it  is  a  slip  surface  across  which  the  pressure 
and  the  normal  velocity  component  n'q  are  continuous  while 
the  density  and  the  tangential  velocity  component  can  have 
arbitrary  jumps.   In  the  particular  case  both  m 
vanishes  and  nxq  is  continuous,  the  slip  surface  is  called 


17 


a  contact  discontinuity  where  only  density  is  discontinuous 
and  there  is  no  relative  motion. 

If  there  is  a  region  of  supersonic  flow  in  the  flow 
field   it  is  well  known  [30]  that  shock  waves  will 
generally  appear.   The  entropy  condition  to  pick  the  right 
weak  solution  is  that  n-q  decreases  across  the  shock. 

For  the  inviscid  flow  model  the  boundary  condition  at 
the  body  is  reduced  to  the  kinematic  condition  requiring 
the  body  to  be  impenetrable   to  the  flow.  Namely,  the  flow 
remains  tangent  to  the  body  surface.   In  mathematical  terms 
it  is  subject  to  the  condition 


H^  =  F   +  q-VF  =  0   on  the  body  surface   F(x,y,t)  =  0 


At  the  trailing  edge  it  requires  that  the  pressure  and 
the  flow  direction  be  continuous.   To  be  specific,  the 
rate  of  change  of  circulation  r,  measured  counterclockwise, 
is  given  by 


^-A^^^^ 


i^    •    ^r.m 


=  -  (t  ^  +  i  ^ 


dt         7   2 
2 
P   ■  T   2~ 
2 


,  >o:r 


dq' 


,  2,  ^qu-^q^^  , 

=  [q  ]    = 2 ^^u  -  ^£^ 
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where  the  velocities  q  and  q^  are  the  upper  and  lower 
velocities  at  the  trailing  edge.   Consequently,  if  the 
circulation  is  to  change  with  time,  neither  the  average 
velocity  nor  the  jump  velocity  can  be  zero.   The  vortex 
sheet,   comprised  of   vortex   filaments,   trailing 
downstream  of   the   airfoil,   is  viewed   as   a 
slip    surface.     In   reality,   the   vortex   sheet 
is  convected  with  the  motion  of  the  fluid  and  rolls  up 
on  itself  due  to  its  self-induced  velocities.   A  consistent 
model  accounting  for  the  roll-up  of  the  sheet  would  add 
greatly  to  the  difficulty  of  constructing  a  boundary 
conforming  coordinate  system.   If  the  convection  and  roll-up 
of  the  sheet  are  ignored,  the  vortex  sheet  may  be  assumed 
to  be  along  the  streamwise  coordinate  surface  that  leaves 
the  airfoil  trailing  edge  smoothly.   The  constraints  applied 
on  it  are  that  the  pressure  and  the  normal  velocity  component 
be  continuous  across  the  vortex  sheet. 

The  appropriate  radiation  boundary  conditions  at  the 
artificial  boundaries  of  computational  domain  are  again 
needed. 

We  remark  that  in  steady  flow  calculation,  the  energy 

equation  can  be  replaced  by  Bernoulli's  equation  for  constant 

2   2 
total  enthalphy  H  =  -^  -  +   ^  „^   =  const.   thereby  reducing 

the  number  of  dependent  variables  from  four  {p,u,v,e)  to 

three  (p,u,v) . 
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1 . 3   Potential  Flow  Equation 

Assuming  that  the  flow  can  be  described  by  a  velocity 
potential,  the  Euler  equations  can  be  reduced  to  a  single 
quasilinear  equation.   This  implies  that  the  flow  is  irrota- 
tional  and  hence  in  view  of  Crocco's  relation  that  there  are 
no  entropy  changes  in  the  flow.   The  entropy  produced  by  a 
shock  is  proportional  to  the  third  order  of  the  shock  strength 
[12].   We  may  assume  that  the  entropy  is  conserved  across  the 
shock  if  we  just  consider  weak  shocks,  such  as  occur  on  the 
surface  of  a  well  designed  airfoil.   This  approximate  model 
should  not  be  a  source  of  serious  error  if  the  Mach  number  of 
the  normal  component  of  the  flow  ahead  of  the  shock  is  less 
than  1.2. 

Let  $  be  the  velocity  potential  with  q  =  V$   the 
velocity  vector.   The  equation  of  motion   Dq/Dt  =  -Vp/p 
leads  to 

$^  +  ^  +  j  ^  =  f (t)  +  constant. 

If  ())  =  $  -  I  f(t)  dt  then  V^  =  V$  and  4).  =  $.  -  f(t). 
We  therefore  call  cj)  velocity  potential  as  well  and 
we  have  the  Bernoulli   equation   for   (J) 

.2    r 

+ 
J 


(1)  *t  ^  V  + 


—  =  constant. 

P 


The  conservation  of  mass  is 
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(2)  P^  +  (P*x^x  "^  (P'^y)y  =  0 

We  take  the  equation  of  state  to  be 

(3)         '    P  =  (^)   P^ 
with    p   =1   and   a   =  1  . 

■^  CO  CO 

In  the  smooth  region  of  the  flow  we  may  eliminate  p  from 
the  above  equation  and  get  a  quasilinear  equation  for  (p. 
The  equation  of  continuity  yields 


1  Dp 
P  DtJ 


=  Vq  =  V'V(()  =  A(}) 
The  Bernoulli  equation,  after  differentiation,  leads  to 


D       -2 


Dt  '"^t   2 


Y-1 
dpi  ^  _  _D_  fP_[ ^  ^  _  ^Y-2  Dp 

Dt 


'*t-V'  --l:   f  --m(^)'-- 


a   Dp     2.  , 
.   =  -  T  Dt  =  ^  ^'^ 

Finally,  we  combine  them  and  get  an  equation 

(4)  ^  ^^t  "^  2  '^^^  =  ^^^* 
or 

(5)  't'tt  +  2u<}.^^  +  2vct)y^  =  (a2-u2)(t.^^  -  2uv't>^^   +    (a^-v^)<i>^^ 

The  shock  conditions  which  will  be  applied  in  the  model 
(1)  ,  (2)   and  (3)   are 

a.  nxq  is  continuous  across  the  shock  which  implies 
())  is  continuous 

b.  (n'q  -  s)'p   is  continuous   which  says  mass  is  conserved 
across  the  shock,  where  s   is  the  shock  speed 
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c.  n-q  decreases  across  the  shock.  This  is  the  entropy 

condition. 
Here   n  is  the  normal  to  the  shock  surface. 

According  to  these  conditions,  a  normal  shock  is  to  be 
modeled   as  a  jump  between  equal  points  of  an  isentropic 
stream   tube.   The  corresponding  change  in  normal  momentum 
is  balanced  by  a  force  on  the  discontinuity.  The  combined 
force  on  the  body  and  the  discontinuity  is  zero  so 
that  the  integral  of  the  pressure  over  the  body  surface 
yields  a  drag  which  is  an  approximation  to  the  wave  drag. 

The  surface  condition  requires  that 

V(})«n  =  v_*n   on  the  body  surface. 

Here   n  is  the  normal  to  the  body  surface  and  v   is  the 
body  velocity  relative  to  the  absolute  reference   frame. 
The  trailing  edge  and  wake  condition  are  basically   the 
same  as  for  the  Euler  equations.   Specifically,  the  rate 
of  change  of  circulation  V   of  airfoil  is  given  by 


dr    d 


d 


(^)        It  =  ai  t  ^  ^^  ^  5t  f^^TE  =  f^t^TE         ■  ^' 

With  the  continuity  of  pressure  and  normal  velocity  component 
across  the  wake,   the  Bernoulli  equation  gives 

(7)        [(J)^]  +  [^]  =  [^^]    +   4>s[<t>s^  =  ° 

where   *   stands  for  the  average  velocity  at  any  point  in  the 
s 

wake.    Thus   the   circulation   can   change   only   if 
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there  is  a  velocity  jump   at  the  trailing  edge.   Hence  a 
vortex  sheet  is  shed  and  the  wake  condition  (7)   expresses 
the  equation  for  the  transport  of  vorticity  downstream. 
We  will  discuss   radiation   boundary   conditions   for 
equation  (5)   in  a  later  section. 


1. 4   Small  Disturbance  Equation 

For     small  disturbance  transonic  flows,  the  flow 
equation  can  be  further  simplified  by  a  perturbation 
method  [1]  .   Namely,  assume    that  the  thickness   to 

chord  ratio   t   of  the  airfoil   under  consideration 

2/3       2 
is  small  in  the  sense  of    x     '^>l-M  <<1,  where  M   is 

the  free  stream  Mach  number.   If  we  expand  the  potential 

cj)  to  the  potential  flow  equation  in  the  powers  of  x  and 

retain  the  lowest  approximation,  we  obtain  the  small 

disturbance  equation 


^l^tt  -^  ^^iKt  =   ^c*xx  -^  * 


where 


and 


yy 


S,  =  M^(k:2/x^/^)  ,   S_  =  M^(k/x^/'^) 

J.         oo  £.  '" 


V   =  (l-M^)/x^'^^-  (y+1)M^(|)  -  (Y-l)M^K(j). 

Q  ^  CO'  '       °°X  oo~-^ 

The  reduced  frequency  k  =  toc/q   is  a  measure  of  the 
degree  of   unsteadiness  of  the  flow  field  since  it  is  the 
ratio  of  the  time  scale  of  the  airfoil  flight  speed  c/q 
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and  that  of  the  unsteady  motion  1/w,  where  c  is  the  chord 
of  airfoil,  w  is  the  frequency  of  the  unsteady  motion  and 
q   is  the  flight  speed.   The  flow  velocity  is  the  sum  of 
the  free  stream  velocity  q^  and  the  gradient  of  ({).   We 

remark   that  (^ ,    t,  y  and  x   have  been  scaled  by 

2/3      w     /  1/3    ,  ^-   n 

CT    q   ,  l/w,  c/T     and  c  respectively. 

The  primary  merit  of  this  approximation  is  that  the 

surface  condition  is  very  simple.   The  surface  of  the 

airfoil  is  transferred  to  the  slit  y=0,  0<x<l, 

which  is  the  mean  surface  approximation  to  the  airfoil 

in  the  new  scaled  coordinate  system.   If  h(x,t)  is  the 

unsteady  displacement  of  the  airfoil  surface  from  the  true 

mean  contour  f(x),  then  the  surface  condition  is 

(b      =    f      +   h      +   h^      on   the    slit      y=0,       0<x<l. 
y  X  X  t  J  ' 

The  wake  condition  is  that  the  jump  of  the  pressure  coefficient 

across   the  wake   y  =  0  ,  1  <  x,  must  vanish  Namely, 

[c^]  =  0  where  c^  =  -  2  x^^-^  (ct)^ +<}>.)  . 
P  P  X    t 
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1. 5  Low  Frequency  Small  Disturbance  Equation 

2     2/3 

For  low  frequency  <  ^   1-m  '\.  t  '   <<  1,  it  is  well 

known  [1]  that  the  small  disturbance  equation  reduces  to 

2^xt    c^xx    ^yy 
where   v  =  (1-M  )/x^        -    (y+1)M^4' 

Q        *        OO''  '  °°X 

The  surface  boundary   condition   and   the  wake   condi- 
tion  can  be  either   that  of  the  small  disturbance 
equation  or  as  follows: 

a.  4)      =f      +h        on      y=0,       0<x<l 

y  X  X  ■' 

b.  [c    ]    =    0  on      y   =    0    ,      1    <   x      where   c      =    -2t    ' -" i;,    . 

p  II  p  ^x 

1.6  When  to  Use  Which  Model 

Each  model  has  its  own  limitations  based  on  the  assump- 
tions used  in  developing  the  flow  equations.  For  example, 
the  low  frequency  small  disturbance  equation  does  not 
describe  high  frequency  motion  well,  the  small  disturbance 
equation  does  not  describe  the  blunt  leading  edge  airfoils 
well,   the  potential  flow  equation   does  not  describe  the 
strong  shock  wave  well,  the  Euler  equation  does  not  describe 
separated  flow  well.  We  briefly  remark  that  when  strong  shocks 
lead  to  separation,  viscous  effects  cannot  be  neglected.  Either 
boundary  layer  correction  equations  or  the  Navier-Stokes  equa- 
tion have  to  be  employed  [10] .  The  consideration  of  turbulence 
is  probably  needed  to  resolve  the  complicated  flow  phenomena 
such  as  buffet  separation,  reattachement,  and  so  on.   Here  we 
will  consider  flows  with  relatively  weak  shocks  which  can  be 
adequately  described  by  the  potential  flow  model. 
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2 .    Basic  Numerical  Concepts 

The  numerical  problem   is  to  find  an  approximate 
solution  accurate  to  within  some  tolerance.   The  most 
basic  and  widely  used  method  to  solve  time  dependent 
problems  is  the  finite  difference  method.   In  this  sec- 
tion we  review  some  basic  numerical  concepts  about  the 
finite  difference  method  and  propose  a  finite  difference 
strategy  with  a  splitting  technique  which  result   in 
unconditionally  stable  schemes  for  the  heat  equation, 
linear  advection  equation  ,        and  wave  equation,   respectively. 
And  we  will   apply   those   ideas   to   construct   an  ADI 
type  scheme  for  the  potential  flow  equation  in  quasilinear 
form  in  Section  IV. 


"-■'  4  -i 
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2.1   Mesh  Spacing 

In  the  finite  difference  method,  one  performs  all 
calculations  on  the  grid  points  of  a  computational  domain 
which  is  of  finite  extent.   Once  the  grid  points  are 
given,  the  resolution  of  the  physical  phenomena   is 
naturally  limited  by  the  mesh  spacing.   To  be  specific, 
we  introduce  some  terms  through  the  definition  [36]  of 
a  Fourier  mode; 

y  =  a  e     ^ 

^  ^  ^i^(x+(co/^)t) 

i  ^(x+ct) 

(1)  =  a  e   ^ 

where   a  is  called  the  amplitude;  o),  the  phase  rate; 
E, ,    the  wave  number,  c  =  oo/^,  the  wave  speed;  A  =  2TT/Cf  the 
wave  length;  ut  +  E,x,    the  phase  angle;  f  =  a)/2Tr,  the 
frequency;  t  =  1/f,  the  period. 

Suppose  we  express  a  function  u(x)  as  a  Fourier  series 

i^.x 

(2)  u(x)  =   I  a.  e   -> 

— oo    -* 

On  a  mesh  system  containing  I  equally  space  points  of 

spacing  Ax,  the  Fourier  mode  of  shortest  wave  length 

resolvable  in  the  system  is   A  .   =   2Ax  ;  the  longest 

■^  mm 

wave   length   is      A  =    (I-l)Ax   =   L.      The   corresponding 

^  max 

wave   numbers   are    E  =   tt/Ax   and   ^    .      =   2it/L. 

^max  min 

So   the  total   number  of  wave  models   resolved  by  this 
mesh  system  is  N  =  (I-l)/2   and  the  part  of  u  which  can  be 
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resolved  by  this  system  is  the  partial  sum 


N      i^  .X 
u  =   y  a  .  e   -'    with   ^  .  =  ^-j— 
_N   ^  IJ    NAx 


In  viscous  flow,  the  diffusion  and  the  advection  for 

iEx 
a  Fourier  mode  u  =  a  e     lead  to 


d^u        p2 

y  — 9  =  -  y  ^  u 

dx^ 


and 


dU  t   ■ r\ 

pu  ^  =  pudOu 


Their  ratio  is  (pu)/(y^).   As  ^  increases  the  diffusion 
becomes  stronger  and  dominates  eventually.   The  mesh 
spacing  should  be  fine  enough  to  understand  the  dissipa- 
tion mechanism.   On  the  other  hand,  for  computational 
efficiency   the  number  of  mesh  points  must  be  kept  to 
the  minimum  required  to  resolve  all  the  significant 
phenomena.   Hence,  in  practice  [32],  a  typical  computational 
domain  consists  of  a  fine  mesh  region  where  viscous  effects 
are  important  and  a  coarse  mesh  region  where  the  flow  is 
essentially  inviscid.   Some  techniques,  for  instance, 
coordinate  stretching  and/or  coordinate  transformations  are 
useful  [14,19,24,41].   Automatic  mesh  system  generation 
techniques  for  flow  about  multiple  bodies  in  a  plane  have 

been  developed  [42,43]. 

V-  '■•-■  ■•  ^■■^■i>('' 
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2. 2   Time  Step   and  Approximation  Factorization  Technique 

Explicit  finite  difference  methods  have  demonstrated 
their  ability  to  solve  a  wide  range  of  flow  problems.  How- 
ever the  size  of  a  time  step  that  a  solution  can  be  advanced 
during  each  step  of  calculation  is  restricted  by  the 
Courant-Friedrichs-Lewy  condition  (CFL  condition) .   The  CFL 
condition  imposed  on  the  time  step  is 

A 
At  <  


lq| 

where  A  is  the  grid  mesh  spacing  and  |q|  is  the  fastest 
propagation  speed  anywhere  on  the  mesh  system.  Therefore, 
the  solution  requires  long  and  expensive  computation  time. 

Unlike  the  explicit  method,  implicit  methods   can  be 
theoretically  stable  for  all  time  step  sizes.   Unfortunately, 
an  implicit  method  in  two  or  higher  space  dimensions 
requires  a  set  of  equations  to  be  solved  at  the  advanced 
time  level   which  is  not  always  easy  to  accomplish  directly. 
Accordingly,  the  splitting  technique  is  introduced  to  yield 
feasible  computational  processes.   We  illustrate  the  split- 
ting technique  on  the  heat  equation  in  two  space  dimensions. 


(1)  ^^  =  <J>    +  <}) 

^t   ^xx   ^yy 

The  finite  differencing  strategy  is  replacing  the 
differential  operator  in  time  D  by  the  forward  difference, 
the  (j)  in  the  right-hand  side  by  the  average  of  0    and  4) 
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and  the  differential  operators  D   and  D   by  the  second 

XX      YY 

order  center  difference  operators  in  x  and  y  respectively, 


Namely, 
(2) 


,n+l  ,n 

(p  .  .  -(p  .  . 

At 


E  -2I+E 

X        X 


Ax' 


J 


^   -^,„— 1    ,n+l , , n 
E  -2I+E   -i  >,  /-(^  .  .  +4)  . 


Ay' 


ID 


ID 


n 


,n 


,n    ^n 


where  E^(},.  .  =   <l>  ^^^  ^  .    ,    similarly  E^^^  =   ^i^j+i- 

The  accuracy  of   the  finite  difference  equations  is  of 

second  order  in  time  and  space. 

We  may  write  the  finite  difference  equation  in  terms 

of  6    =  E  -2I+E~''",   5    =  E  -2I+E~"''  ,  with  p  =  At/2Ax 
XX    X     X     yy    y     y 

2 
and  q  =  At/2Ay   as  the  equation 


(3)  fl-p6   -q6  ]^ 
^   "^  XX  ^  yy*  ^y 


n+1 


fl+pS   +q5   ](!)" 
>-  ^  w  ^  yy  y 


XX 


The  idea  behind  the  splitting  technique  is  to  generate 
a  perturbation  of  the  above  equation  that  permits  a  simpler 
computational  process.   Namely,  we  may  factor  equation  (3) 
as  follows. 


(4)   (l-p5_)(l-q6yy)(t^+l   =   (l4-p6^J[l+q6^y)c^J 


xx^ 


Here,  we  add  a  term   (At  /4)<}'j^   ^^  of   third   order  in  time 
to  the  equation  (4) .  The  von  Neumann  stability  analysis  shows 
that  the  scheme  is  unconditionally  stable  which  means  there 

is  no  restriction  on  the  time  step  At  to  the  spacial  steps 

'^k   i  (niX"i~n,v)  • 
Ax  and  Ay.   Indeed,  substituting  (f)  =  (J)   e      ^      into 
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equation  (3)  we  obtain 


*l  - 


1  -  2p(l-cos  K) 
1  +  2p(l-cos  i) 


1  -  2q (1-cos  n) 
1  +  2q(l-cos  n) 


with  E,   =  m   Ax      and  n  =  n  Ay.    By  the  fact  that  both  p 
and  q   are  positive  we  conclude  that  the  right  hand  side 
is  less  than  1.   This  shows  the  amplification   1 4>  I  is 
bounded  by  unity  without  any  restriction  on  p  and  q. 

The  algorithm  for  the  solution  of  equation  (4)  con- 
sists of  three  easy  stepsi 

X  =  (l+q5   )())" 

yy  y 

(l-p5   )Y  =  (l+p6   )X 

'^  XX  '^  XX 

(l-q6   )(f)""^^  =  Y 

Each  of  the  last  two  steps  requires  a  3-diagonal  matrix 
solver  which  is  not  expensive  at  all  and  can  be  found  in 
any  standard  numerical  method  book  [23,41]. 

It  is  worthwhile  noting  that  equation  (4)  can  be  taken 
to  represent  an  iterative  procedure  which  converges  if 

d) .  .  =    (h .  .    =   (h  .  .      for  sufficiently  large   n. 

Then,  equation  (2)  is  reduced  to  the  standard  five-point 
difference  approximation  of  the  Laplace  equations.  In  this 
case  the  quantity  At  can  be  viewed  as  an  iteration  parameter 
and  may  be  varied  from  iteration  to  iteration  to  optimize 
the  convergence  of  the  process. 
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Equation  (4)  can  be  rewritten  as 

(1-  ^  D   )(1-^D   )((^"-rl-<},".)   =   At(D   +D   )({.'?. 

2   XX     2   yyiD   ID        XX  yyiD 
or 

(5)    (a-D   )  (a-D   )( (J)?"^-"- -<])?. )   =   2a(D   +D   )())?. 
^■^'  XX     yy  ^1]   ^1]  XX  yy   iD 


It  falls  in  the  following  general  form  [20,25], 


(6)  Nc"  +  a3R'^  =  0 


which  is  used  to  solve  the  steady  differential  equation  Ij)  =  0 . 
Here,  c   =  (})     -  cj)   is  the  correction,  R  =  L(})   is  the 
residual  which  measures  how  well  the  finite  difference  equa- 
tion is  satisfied  by  the  nth   level  solution  c()  ,  u  is  a 
relaxation  parameter  and  N  is  chosen  as  a  product  of  two 
or  more  factors  indicated  by 


The  factors  N,  and  N„  are  chosen  so  that  (1)  their 
product  is  an  approximation  to  L,  (2)  only  simple  matrix 
solvers  are  required,  and  (3)  the  overall  scheme  is  stable. 

This  type  of  implicit  scheme  has  been  found  very  power- 
ful in  the  calculation  of  steady  flow.   We  remark  that  the 
parameter  a  in  the  equation  (5)  can  be  replaced  by  some 
lower  order  differential  operator  to  speed  up  the  convergence 
rate  as  well  as  to  introduce  damping  which  is  needed  in  the 
multigrid  technique  [26] . 
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2 .  3   Artificial  Dissipation  and  Upwind  Dif ferencincr 

In  inviscid  flow  calculation,   a  scheme   that   seems 
stable  for  shock  free   flows  sometimes  blows  up  when  it 
is  employed  to  calculate  shock  waves.   This  is  due  to  the 
fact  that  using  some  difference  formulas  across  a  disconti- 
nuity can  lead  to  oscillations  which  may  grow.   To  remedy 
this,  the  well  known  shock  capturing  technique  is  to  add 
to  the  inviscid  flow  equation  a  proper  amount  of  artificial 
dissipation  to  simulate  the  physical  dissipation  in  the 
shock  layer  and  to  provide  the  necessary  damping  for  large 
wave  number  disturbances   so  that  the  shock  wave  is  smeared 
out  over  several  mesh  points  [28] .   Namely,  if  we  model  the 
physical  solution  by  the  inviscid  flow  equation 

(1)  u^  +  V'f(u)  =  0 

For  shock  calculations,  we  look  at  the  solution   u  of  (1) 
as    limit   of  the  viscous  flow  equation 

(2)  u^  +  V«f(u)  =  V(eVu) 

where   e  is  positive  and  is  of  the  order  of  the  mesh  spacing. 

Equation  (2)   is  of  diffusion  type  and  the  solution 
can  be  shown  to  exist  [31].   Suppose  that  the  solutions  u(e) 
of  (2)  tend  to  a  limit  u  boundedly  almost  everywhere  as  e  ->■  0. 
Then,  u^(e)  tends  to  u^  ,  Vf(u(e))  to  Vf (u)  and  V{e  V'u)  to  0 
in  the  distribution  sense.   This  says  that  u  satisfies  (1) 
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in  the  distribution  sense  which  is  equivalent  to  saying 
that  u  satisfies  the  conservation  law  in  the  integral  form. 

The  artificial  viscosity  can  be  viewed  as  a  kind  of 
truncation  error  exhibited  by  the  approximation  to  the 
differential  equations.   It  may  be  either  in  explicit  or 
in  implicit  form.   We  consider  the  artificial  viscosity 
introduced  by  upwind  difference   for   the   advection 
equation 


(3)  ^^  +  u(|)^  =  0 

The  finite  difference  approximation  for  the  case  u  >  0  is 

, n+l_ , n 

/  4  \    i     i  ,   u   /,,n+l  ,n+l,  ,  ,,n  ,n   ,1      _ 

(4)   Al^  +  2A^   (<J>i   -*i_i)  +  (<J>i-<t>i.i)    =   0 


The  von  Neumann  local  stability  analysis  is  to  substitute 

^k   imx 
a  Fourier  mode  <})  =  (})   e     into  equation  (4)  .   This  leads 

to 

1  -  p(l-cos  C)  -  ip  sin  C 

<  1 


1 +  p(l-cos  C)  +  ip  sin  ^ 


with  p  =  uAt/2Ax  >  0  and  ^  =  mAx.   It  is  trivial  that  the 
scheme  is  unconditionally  stable. 

In  equation  (4) ,  we  did  add  an  artificial  viscosity 
implicitly  through  the  upwind   differencing  in  x .   We  can  see 
it   explicitly  by  Taylor  series   expansion.    Equa- 
tion (4)  is  equivalent  to  the  equation 

(5)        <t>^   +   u(^^   =   u   ^   4)^^  +  0(At^,Ax^)  . 
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The  extra  term  u(Ax/2)G    is  the  leading  term  in  the  trunca- 
tion  error  and  is  referred  to  as  the  artificial  viscosity. 
To  discuss     the  diffusion  and  dispersion  properties 
of  this  finite  differencing,   we  first  derive  the  dispersion 
relation  of  the  differential  equation  (3) .   Substituting 

,       i(wt+^x)   •   .     .V.  .  •      r-,\  •   1  J    .  v.       1   i.  • 

0  =  e         into  the  equation  (3) ,  yields  the  relation 
(jj+u^  =  0.    CO  is  a  real  number   so  that  there  is  no  damping 
of  any  wave  mode   and  all  waves  have  the  same  phase  speed  u. 

Next,  we  apply  the  same  Fourier  mode  to  the  viscous 
differential  equation 

(6)  *^  +  u*^  =  lu|  ^<J>^^  . 

It  has  the  dispersion  relation 

to  +  u?  =  |u|  ^  C  1 

So  a  solution  of  equation  (6)  is 

i^(x-ut)    -t|u|^  ?^]t 
())   =  e         '6 

The  magnitude  of  the  damping   increases   with  the  wave 
number   ^   and  the  velocity  u.   Hence,  the  effect  of  arti- 
ficial viscosity  is  to  introduce  larger  dissipation  for  both 
the  larger  wave  number  mode  and  the  faster  flow  region. 
We  see  that  there  is  no  dispersion  up  to  the  first  order 
approximation.   However,  if  we  add  an  extra  term  of  (()     to 
the  right  of  equation  (6)   then  dispersion  occurs.   This 
means  that  different  frequency  waves  propagate  with  different 
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speeds  in  the  flow  field. 

The  upwind  differencing  has  played  a  very  important 
role  in  transonic  flow  calculations.   The  main  purpose  is 
to  exclude  the  expansion  shock. 


2 . 4   Conservative  Finite  Difference  Schemes 

The  main  idea  behind  the  use  of  conservation  form  is 
the  fact  that  if  the  difference  equation   to  the  differential 
equation  in  conservation  form  is  again  in  conservation 
form,  the  solution  of  the  finite  difference  equation  satisfies 
the  Rankine  Hugoniot  jump  conditions  automatically  [30,39]. 

The  differential  conservation  form 


(1) 


u   +  div  f(u)  =  0 


can  be  derived  from  the  more  general  integral  form 


Ih 


dx 


t    t 

+ 
s    s  9R 


J 


f-n  ds  dt  - 


(2) 


s  R 

t 

s  R 


u   dx  dt  + 


u   dx  dt  + 


11... 

s    9R 


ds  dt 


t 

r 


V-f  dx  dt 


s  R 


=  0 


which  says  that  the  change  in  the  amount  of  a  substance  with 
density  u  contained  in  the  region  R  of  space  under  considera- 
tion  is  due  to  the  flux   f   of  that  substance  across  the 
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boundary  3R   from  time  s  to  time  t. 

The  conservative  finite  difference  approximation  is 
then  defined  having  the  form 

^n+l 


(3) 


I  1 


, n+1   n-1 

3r    gr 


2At 


1    F 


n-1 


8R 


dt  =  0 


which  simulates  the  integral  conservation  form. 

Our  differencing  strategy  for  the  flow  equation  in 
conservation  form  yields  the  finite  difference  equation 


(4)   I    I 


n+1    n-1 

rUr,       -    Uo 

Br     Br 


^R^ 


2At 


(  1    Fj"-^'+(  I    F^)"-' 
9R  "        8R  "" 


n+1 


=  0 


The  question  is  how  to  solve  for  u    for  this  large 
nonlinear  system.  Some  linearization  for  F   is  needed 
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III.   POTENTIAL  FLOW  EQUATION 

In  the  steady  inviscid  transonic  flow  calculation, 
the  nonconservative  form  method  agrees  well  with  wind 
tunnel  pressure  data  all  the  way  up  to  the  onset  of 
buffet  [18] .   On  the  other  hand,  for  the  conservative 
form  method,  the  agreement  is  less  satisfactory  and  the 
adequate  correlation  with  experimental  data  seems  to  be 
achieved  by  making  correction  with  boundary  layer  shock 
wave  interaction.   For  mesh  sizes  of  practical  interest, 
instead  of  doing  a  better  simulation  by  combining  a  finer 
scale  model  of  boundary  layer  shock  wave  interaction  with 
conservative  transonic  equations,  we  pick  up  the  noncon- 
servative quasilinear  potential  flow  equation  as  our  model 
and  develop  a  computer  code  for  it. 

We  first  discuss  the  characteristic  surfaces  of  the 
equation  and  explain  the  domain  of  dependence  for  super- 
sonic points.  Then,  we  give  a  set  of  radiation  boundary 
conditions  which  is  shown  to  be  very  satisfactory  with 
the  numerical  scheme  we  propose  in  Section  IV.  And, 
finally,  we  introduce  the  coordinate  transformation  such 
that  the  airfoil  is  fixed  on  a  portion  of  coordinate  line. 
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1.    Characteristic  Surface 

It  is  helpful  to  know  the  characteristic  surface  of  the 
flow  equation  on  which  the  wave  front  along  with  informa- 
tion  is  propagated  throughout  the  flow  field.   Let  s  and  N 
be  coordinates   in  the  local  stream  and  normal  directions 
respectively.   The  direction  cosines  of  s  are  u/q  and  v/q. 

d)    and  *,,.,  can  be  expressed  locally  in  terms  of  the  actual 
ss      ^NN 

coordinates  as 

12  2 

^  =    — t"  (u  (})    +  2UV(})    +  v  (()   ) 

^ss    2    ^xx      ^xy      yy 

12  2 

4'»T>T  =  —^    (v  (b    -  2uv())    +  u  d)   ) 
^NN    2    ^xx      ^xy      yy 

The  potential  flow  equation  in  Cartesian  coordinates  locally 
aligned  with  the  natural  coordinate   system  (s,N)  can 
be  written  as 

*tt  +  2q*3t  =  (^^-^^)*ss  +  ^^*NN  • 

The  characteristic  surface   satisfies  the  equation 

2   2 
{q^-a^)t^  -  2qst  +  s^  -  (^  "^  ) N^  -   0. 

a 

As  shown  in  Figure  2,  on  the  (N,t)  plane,  the  character- 
istic equation  is  reduced  to 

a^t^  -  N^  =  0     or     (N-at) (N+at)  =  0  . 
The  disturbance  propagation  speed  is   a. 

On  the  (s,t)  plane,  the  characteristic  eauation  is  reduced 

to 

,  ^    n2     2.2 

(qt  -  s)   =  a  t 

or 

(s-(q+a) t) (s-(q-a) t)  =  0  . 
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The  particle   speed  is  q,  the  upstream  propagation  speed 
is  q-a   and  the  downstream  propagation  speed  is  q+a .  Thus 
the  disturbance  information  is  propagated  by  the  Doppler 
shifted   sound  wave  velocity.   For  transonic  flow,  the 
particle  and  downstream  waves  quickly  travel  away  from 
the  airfoil   but  upstream  waves  remain  in  the  vicinity 
of  the  airfoil  for  a  much  longer  time.   The  slow  waves 
force  a  slow  approach  to  a  steady  state  solution,  while 
the  fast  waves  stipulate   a  small  time  step  by  the  CFL 

condition  At  <_   Ax/ (q+a)  . 

2   2 
If  a  new  time  coordinate  T  =  t  +  qs/(a  -q  )  is 

introduced,  then  the  potential  flow  equation  can  be 

expressed  as 

2 

^^^-"i^^^ss    -^  ^^*NN  =  ,    2       2/  *TT 

(a  -q  ) 

2     2 

So  for  supersonic  points,  a   <^  q   ,  T  is  a  space-like 

direction   and   s   is  a   time-like  direction.   This  means 
the  differencing  in  the  s-direction  should  be  retarded  in 

the  supersonic  region   in  order  to  have  the  right  domain  of 

2     2 
dependence.   For  subsonic  points,  a   >  q   ,  s  is  a  space- 
like direction  and  T  is  actually  a  time-like  direction. 
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(a) 


(b)   subsonic  case,  q  <  a 


/  1 

/  q-a 


q+a 


(c)   supersonic  case,  q  >  a 


Figure  2.   Characteristic  Surface  of  the  Potential 
Flow  Equation  in  Ouasilinear  Form; 
(a)  (N,t)  plane;  (b) , (c)  (S,t)  plane. 
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2 .    Computational  Boundary  Conditions 

Problems  of  transonic  flow  field   are  usually   posed 
in  the  exterior  of  the  body  which  is  an  unbounded  domain 
[4].   Owing  to  the  finite  storage  capability  of  the 
computer,  the  numerical  computations  require  that  the 
computational  domain  be  finite.  The  proper  boundary  condi- 
tions must  be  developed  at  these  computational  boundaries 
so  that  the  computed  solution  closely  approximates  the 
free  space  solution  which  exists  in  the  absence  of  these 
computational  boundaries  [15,16,17]. 

For  steady  state  calculations  in  transonic  flow, 
coordinate  mapping  techniques  are  a  traditional  and 
effective  way  of  handling  these  computational  boundary 
problems.   The  reason  for  the  success  of  coordinate  mapping 
techniques  lies  in  the  fact  that  the  steady  state  far  field 
asymptotic  behavior  is  given  by  a  regular  algebraic 
singularity  without  oscillation.  For  genuinely  unsteady 
transonic  phenomena,  the  solution  of  flow  equations  usually 
possesses  a  strongly  oscillatory   transient  behavior  and 
the  far-field  asymptotic    behavior  is  an  oscillatory 
singularity.   The  standard  coordinate  mapping  technique  is 
not  adequate  to  resolve  this  problem.   It  must  be  supplemented 
by  a  set  of  proper  boundary  conditions  at  the  computational 
boundaries.   "  ■   '-     ,"   "    i   ^.  ;• 

In  this  section  we  will  give  a  set  of  radiation  boundary 
conditions  for  the  potential  equation  in  the  Cartesian 
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coordinate  system  and  in  a   later   section  we  will  give 
its  corresponding   form  in  the  computational  domain. 
In  the  physical  domain,  the  computational  region  for  an 
airfoil  in  two  dimensions  is  depicted  as 

R- 


Q 


2 


Figure  3.   The  Typical  Computational  Region  for  an  Airfoil. 

The  design  of  effective   far  field   radiation  boundary 
condition  depends  on  the  wave  propagation  properties  of 
the  flow  equation.   We  consider  the  potential  flow  equation 


(1)  <J)^^+  2u4)^^+  2v<^y^=  (a^-u2)({)^^-  2uv(J>^y+  (a^-v^).})^^ 

For  a  plane  wave  $  =  e^  ^'^^'''^^"^^^^  to  satisfy  equation  (1)  , 
it  requires  that  its  wave  information  satisfy 

(2)  u?   +  2u£;u  +  2vnto  =  (a^-u^)C^  -  2uv^ti  +  (a^-v^)n^ 

2     2   2     2 

or  (u+  u^  +  vn)  =  a  (^  +  n  ) . 


A  boundary   condition  on   the  upstream  wall,  R, 
boundary,  which  annihilates  the  upstream  propagating 
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wavelet  is  given  by 


(4)  w  +  u?  +  vn=  a?/l+n^/^^ 


Recalling   the   dual  relationship   between   io),  i^ ,    in 
and  D.  ,  D   ,  D    respectively,   the  equation  stands  for 

a  nonlocal  condition.   By  the  first  approximation  of 

12       3 

/1+x  =  1  +  1/2  X  -  —  X   +  0(x  )   we  get  the  first  radiation 

8 


condition  for  R,  boundary,  namely,   03  +  u^  +  v^  =  a^ 
which  leads,  after  Fourier   transformation,  to  the  condi- 
tion 

(5)  4)^  +  (u  -  a)(})^  +  v{j)y  =  0 

Similarly,  we  can  derive  the  artificial  boundary  condi- 
tions for  the  R„  ,  R-  and   R.  boundaries.   At  their  inter- 
section points   P,  ,  P„  ,  Q,  and  Q„  ,  we  use  the  average  of 
the  corresponding  conditions,  and  have  the  following  general 
formula 

(6)  4)^  +  u4)^  +  v(t)^  =  0 

with   u  +  IV  =  (u  +  iv)  +  ae    where  B  =  ~  x  '  J  '  ~4  '  "4 
at  Q2  ,  P2  ni  Q,  and  Q,   and  6  =  0,2'   tt,  — 2  ,  on 


^2  f  Rt  f  R-i  and  R.   respectively. 
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3.    Coordinate  Transformation  Technique 

When  the  body  surface  crosses  the  coordinate  lines 
it  is  difficult  to  satisfy  the  physical  boundary  condi- 
tions.  This  is  particularly  the  case  near  the  leading 
edge  of  the  modern  supercritical  wing  section  where  the 
surface  has  a  high  curvature  and  the  flow  is  sensitive 
to  small  variations  in  the  shape  [41] .   For  the  rigid  body 
motion  of  the  airfoil  the  treatment  is  facilitated  by  the 
use  of  a  moving  sheared  parabolic  coordinate  system  in 
which  the  body  contour  coincides  with  a  segment  of  coordi- 
nate line   and   the  whole   mesh   system   is   moving 
with  the  wing  section  so  that  the  relative  position  of  grid 
points  is  kept. 

We  describe  the  moving  sheared  parabolic  coordinate 
system  as  follows  [3,25]: 

3.1   Coordinate  System 

First,   we  consider  the  physical  plane  to  be  described 
in  a  Cartesian  coordinate  system  (x,y) ,  and  the  airfoil 
attached  coordinate  system  in  Cartesian  coordinate  system 
(x*,y*) .   Let   the  origin   of   {x*,y*)  system  be  at  the 
singular   point  of  the  parabolic  mapping  which  unwraps  the 
airfoil   and  will  be  described  in  the  next  step.  If  the 

■i  (  ■JT  —  O  ) 

flight  velocity  of  the  airfoil  is  M^e  ^      at  time  t,  then  the 


position  of  the  origin  of  the  (x  ,y  )  system  can  be  described  as 

j^^^i(Tr-e)  ^g^   j^  ^^^      angle  of   attack   of 
0 


* 
0  0  = 


J 
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the  airfoil  at  time  t  is  a,  then  the  x*-axis  on  which  the 
airfoil  chord  lies  will  have  an  angle  -  (a  +  9)  with 
respect  to  the  x-axis.   Their  relation  can  be  seen  in 
Figure  4   and   described  as  the  relation 


,-,\        /  .•  %    f  w  /  X   i{iT-9(s))  ,   ./*.•*%   -i(6+a) 
(1)   (x+iy)  =   M^(s)  e  ^    ^  ' '  ds  +  (x  +iy  )  e   ' 


Figure  4.   Frames  of   Reference, 
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Second,  we  unwrap  the  airfoil  by  introducing  the  square 
root  mapping 

(2)  2(x*  +  iy*)  =  (x^  +  iy^)^ 

which  maps  the  entire  airfoil  contour  to  a  shallow  bump 
near  y,  =  0  /  as  shown  in  Figure  5b. 

Third,  if  we  denote  the  height  of  the  bump  as  y,  =  s(x,), 
then  the  shearing  transformation 

(3)  X  +  iY  =  x^  +  i(y^  -  s{x^)) 

reduces  the  airfoil  contour  to  a  portion  of  the  line  Y  =  0. 
Fourth,  we  stretch  the  coordinate  line  by  the  stretch 
mapping  to  render  the  computational  domain  finite.  The 
stretch  mapping,  for  instance, 

bY 

(4)  Y  =  ,   0  <  a  <  1 

(l-Y-^)^        -   - 

will  map  the  infinite  lines  Y  =  +  oo   toY  =  +  l. 

Fifth,  avoiding  discontinuities  at  the  trailing  edge  of  the 
wing  section,  the  branch  cut  is  contined  smoothly  downstream.  In 
physical  space,   the  continuation  is   represented  by 


(5)       y  =  Yte  ■*■  ^f^te-^*l 


where   x  is  the  mean  of  the  upper  and  lower  surface  slopes 

_    _        -* 
at  the  trailing  edge  (x.  ,y   )  and  x   is  a  suitably  chosen 

scaling  constant  (usually  taken  as  the  ordinate  of   the 

local  quarter-chord  point) . 


£n 

[ 

-    -* 
x-x 

xte-x* 

] 

[- 

_    _* 
x-x 

] 

_          _* 
X.      -X 

te 
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A^ 

^ — 

-^ 

B 

^ 

y^ 

C 

(a)  Cartesian 
coordinates 


(b)  Parabolic 
coordinates 


(c)  Sheared 
parabolic 
coordinates 


vv. 


Figure  5, 
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3. 2   Flow  Equation 


The  transformations  (1)  and  (2)  are  conformal.  We  will 
write  the  flow  equation  on  the  (x, ,y, )  coordinate  system, 
and  use  the  chain  rule  to  convert  the  equation  into  the 
(x,y)  system.   Several  key  formulas  are  written  down  for 
reference.   We  begin  with  some  notation. 

Let  z  =  X  +  iy,    z   =  x*  +  iy   , 


z^=   x^+    iy^,   Z   =  X  +  iY 


,   Z  =  X  +  iY 
Then  the  mapping  (1)  and  (2)   may   be   expressed   as 
the  following  compact  forms. 


z  =   I  M^(s)  e 
0 


i(TT-e(s))  ^g  _^  ^*  ^-i(6(t)+a(t)) 


z^    =    2z 


The  modulus  of  the  mapping  function  to  the  z,  plane   can 
be  evaluated  as 


idz 


The  velocity  components  in  the  (x,  ,y-.  )  system: 


u  = 


H 


V  = 


H 


The  chain  rule  gives  the  relation  for  4)   in  the  (x,  ,y,  ) 
system: 


*. 


dx, 
z  fixed     1    ^1 


+    <t> 


dyi 


z  fixed 


dt 


z  fixed 
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wh 


ere  dx,/dt  and  dy,/dt  can  be  found  as  follows.  Since 


2    ^  *    ^  i(e+a)  J      f^  ^,    i(Tr-6) 


z,   =  2z   =  2e  '    M  z  -     M^  e  '    '  ds 

we  take  differentiation  with  respect  to  time  t  and  hold 
z  fixed.  Hence 


' z  fixed  1 


We  remark  that  (j).  =  4)„   -  (v^-V){j)  if  v_,  is  the  relative 

t       i  ^         K  K 

velocity  of  the  origin  of  the  (x,,  y, )  system  to  the  (x,y) 

syatem  [36]  .   Then,  we  conclude  that 

dx,      dy^ 
^R  =  -  («  dF-  '  «  dt-^  • 

The  same  differentiation  applied  twice  on  z^  yields 


d  z,  z,  Zt  rdM 


'1 


dt2 


r  fi(^t+«tt)^  ^  ildt^^  M,(ia^)p^ 


z  fixed 

(Qt^^t^'       ^*   i2a 
4—^1  ";|"     ' 

which  will  be  needed  in  the  evaluation  of  the  following  term: 

2         2 
dx-j^         dy^       d  x^      d  y^ 

dx,  dy^^   dxj^ 

■.    +  ^\t^  -^  *x^x^  dir  -^  \yi  d^)  dt-    .  ;   ■ 

dx^  dy^   dy^ 

+  f*y^T^  +  *x^y^  d^  -^  ^^y^L   dt-J  dt" 
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Similarly,  the  chain  rule  gives 


dx^  dy         dx,  dy 

^x^t    ^^T^    ^1^1  '^^       1^1        ^x^'-dt  -"x^   ^yj^^dt  ^x^ 

dx^         dy^        dy,  dy 

^y^t    ^y^T^    ^x^y^^  dt     ^y^y^  dt    ^y^^  ^dt  ^y^    ^y^  ^dt  ^^ 

Recalling  that 

I  (V(^V)(|V((,|2)  =  J^.  (u^A      +  2uv((,      +  v^A 

z  ^^    \        x^x^       x^y^      y^y^ 

2^  2  1 

,u  +v  ,  ,     ,     . 
-  ( — ij — )  (uxj^  +  vy^)| 


and 


A(J)  =  ^  {(j)      +  (|)     } 
h2   '^'^l^l     ^1^1 


Finally,  we  can  write  down  the  potential  flow  equation 
in  the  (x, ,y, )  frame  as  the  partial  differential  equation 

2 

U  V  U  U  V 

Tj^T^      H   x^T^      H   y^T^     xj^x^  ^^      x^y^   ^^ 
2         2 
Vl  ?    \iZ?"  ^  '  H  ^d^^x^  ^  IT  (dt-)yj 

-'I  ^  dt  x,  y,  ^      H 

2 

=   ^  {(j)      +  (J)     } 
H^   ^1^1     ^1^1 

dx^  dy^ 

where  u   =  u  +  H  j-r—  and  v   =  v  +  H  j^—  . 
r         dt         r         dt 
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The  shearing  an'?   the  stretching   trans fornations 

will  further  bring  the  flow  equation  in  much  more 

complicated  form.   We  will  not  write   down   the  flow 

equation  in  the  (X,Y)  frame  here.  Instead,  we  write  the 
useful  formulas 


())  =  (j)   -  s  A 

X,  X     X  Y 

<^-  =  *„ 

^1  ^ 

2 

d)  =  d)    -  2s  d)    +  s  (j)    -  s   6 

x,x,  XX      X^XY     X^YY     XX  ^Y 

d)  =  d) 

^1^1  ^^ 

*y^T^  =  *YT 

and  ,  J,   dX 

*X  =  *J  dX 

*Y  -  *5  dY 

:   ♦XX  -  ♦xX  @'   *  *J   f^ 

.  ^-  *x.  =  *Sv  ©© 
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3.3   Body  Surface  Condition 

The  velocity  observed  in  the  (x,,y,)  frame  is 
q^  =  ^^r'^r^  =  V(j)  -  V  .   Thus,  the  nonpenetrating  surf 
condition   requires  that  q  •  n  =  0  ,  which  leads  to 


ace 


2f  <^^l   ^^il 


3.4   Wake  Condition 

The  zero  pressure  jump  in  the  wake  which  lies  on  the 
portion  of  the  singular  line  along  which  the  airfoil  is 


opened  up  can  be  [21]  expressed  as 

dx 
T 


f\^  ^  dt^  ^\^    ^   1  ^\^    =  0 


where  u  is  the  average  of  the  upper  and  lower  wake  velocities 


3. 5   Computational  Boundary  Conditions 

The  computational  domain  is  depicted  as 


lower  wake  upper  wake 

airfoil 


O, 


Figure  6. 

The  radiation  boundary  condition  is  of  the  form 


dx. 


^y^ 


dt~  *x^  "^  dt   "^y 


4)   +  u(f)   +  v(i)   =  0 


^1 


where   u  and  v  are  defined  as  before. 
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IV.       NUMERICAL   METHOD 

In    this    section,   we    apply    the    idea   introduced  before 
to   design   our   numerical   solver    for  potential    flow   equation 
in   quasilinear    form.       First,    we    use    type    dependent    dif- 
ferencing  to   introduce    the   proper   amount   of   dissipation 
into    the    finite    difference    approximation    such    that    the 
scheme    is    stable    and   shock   waves    are    captured    auto- 
matically.     Then,    we    factor    the    finite    difference   equation 
into   one-dimensional    factors    so    that  we    are    able    to   solve 
the   equation   efficiently  by   employing   a   5-diagonal   matrix 
solver.       Since    the    disturbance   of   potential    flow   is   prop- 
agated by    the    total   effect  of    advection   and  wave   propaga- 
tion,   we   examine    the    stability   of   our   method   to    two   linear 
models:    advection    and  wave   equations.       The    local   stability 
analysis    shows    that   our    finite    differencing   strategy    and 
approximate    factorization    technique    result   in   uncondition- 
ally   stable    schemes    for    these    two   models. 

1.    Finite   Difference    Scheme 

The    finite    difference    scheme    is    a    time   marching 
alternating   direction   implicit   scheme.      Before  we  write 
down    the    differencing   strategy  we   need   the    following 
convention:  ... 
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(})M    =    (})       -    4* 

D   =    central   difference 

-«- 

D   =    upwind   difference 

D   =   Type    dependent   difference 

To  be   specific,    we    define    the   operators    in   x-direction. 


D    (t).     = 


*i+r*i-i 


2Ax 


4).-<j).     ^ 

r  u*  (  ^    ■^"■)     if  u  >  0 

Ax 


uD    (t).     =   / 
x^i 


4).    ,-(|>. 
l^uM^i—^)      if  u   <   0 


'^  D    for  subsonic  point 
D   =  {  ^^ 

V  D    for  supersonic  point  . 

1. 1  On  interior  point  of  computational  domain 

The  potential  flow  equation  in  Cartesian  coordinates 
locally  aligned  with  the  natural  coordinates  system  assumes 
the  canonical  form  [21] 

with  q  =  (u,v)  the  velocity,  s  and  N  are  coordinates  in  the 
local  stream  and  normal  directions. 
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Basically,  the  velocity  components  u  and  v  at  the 
grid  point  are  evaluated  by  using  central  difference  at 
time  level  n.   The  (J)   term  in  the  Bernoulli  equation  is 
evaluated  by  backward  difference   in  the  natural   way. 

1.11  For  (})    term,  central  difference  will  be  used  for 
temporal   derivative, 

^^  (At)^  (At)^ 

1.12  For  contributions  to  d)  ^,   upwind  differences 

St 

will   be    used    for    all    spatial    derivatives    and   central 

difference  will  be    used    for     temporal   derivagives. 


2g*st  =    2U4,^^   +    2vct,^^ 
=     (2uB   +2v6    )&>^ 

X        y     t 

^n+l    ,n-l 

=     (2uS    +2vB    )  {^—.fri ) 

X  y  2At 

=    -V(uD    +vD    )  ((f,N+(})M). 
A  t         X        y 

1.13   For     contributions    to   ({>       ,    type    dependent   differences 

will   be    used    for   all    terms.       The    term,   (p      is    si±)Stituted 

^^      ^-L.  £      n+1         ,      n-1      .  ,n        it>        +4) 

by    the  mean   of   ({)  and   cf)         ,    i.e.  ,    4)      =   -^ 2^ , 

,  *SS     =     "T^^^^XX-'^^^^^y+V^V^ 

=    ^{u^D      +2uvD      ^v^D   „)  (^ i^ ) 

2  XX  xy  yy  2 

1,2-       . ^       -       ,2-       «   .*N-*M+2*". 

=   — :r(u   D      +2uvD      +v   D      )  [^ ^ ^^^—). 

2 '         XX  xy  yy  2 
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1,14    For     contributions    to   *„„ ,    central    differences  will 
NN 

be    used    for   all    terms.       The    term  <^      is    again    replaced  by 
the   mean    of    cj)  and   cf) 

19  O 

(i)i.i«   =    — :r(v"'(i)      -2uvd)      +U"'*       ) 
^IIN  2         ^xx  ^xv  yy 

q  ^  J  J 

=    —{vD      -2uvD      +u   D       )  (^ -Ti ^^^—)  . 

2         XX  xy         yy  ^ 

Finally,    the    finite    difference    approximation   can  be 
written   as 

^^^  +    A(uO    +vD    )  ((|)N+(j)M) 
(At)2         At        X        y      ^ 

2  XX  xy  yy  2 

^    4    (v^D      -2UVD      -Hu^D       )  (^ImW: 
2  XX  xy  yy  2 

or 

2         9  1 

[l+AtuD^+AtvD^-(^-^)  (^)  (u-5^^+2uvD^+v^5^y) 

-  ^(^)(v'D^^-2uvD^^H-u2D^^^]cfN 

(At)2   [(^)  (u2d^^+2uvD^^    ^v^B^^)    -h    -^(v2D^^-2uvD^y+u^Dyy)]*^ 
~1  4. 

0         9         2 

+  [i-(At)(uD^+vby)  -  i^(^::2a_)  (u^B^^-f^uvB^+v^B^y) 

2       2 

-  (At)    ^a    J  ^^2^       -2uvD      +u^D„)](})M  (1) 


2  2  XX  xy  yy 
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The  discretization  errors  associated  with  the  finite 
difference  approximation  is  of  first  order  in  space  and 
second  order  in  tine.   The  leading  error  terms  in  the 
space  derivative  introduce  the  desired  shock  viscosity. 
The  system  of  algebraic  equation  generated  by  the  equation 
(1)  is  large  and  cannot  be  solved  efficiently.   However, 
this  equation  can  be  factored  within  the  same  order  of 
accuracy  in  time  and  space  by  the  spirit  of  splitting 
technique.   The  following  factorization  has  been  tested 
and  found  to  be  numerically  stable  with  time  steps  much 
larger  than  the  time  step  allowed  by  the  CFL  condition 
for  explicit  methods.   Let  M   =  q  /a 

and 

Then,  the  approximate  factorization  of  the  equation  (1)  can 
be  written  as 

L  'L  AN  =  RHS  (2) 

x   y^ 

This    factorization    reduces    the    large    complicated  matrix 

inversion   problem   to    two   one-dimensional   problems.      The 

algorithm   can  be    expressed    as  ' 

Tl    X    =    RHS      ■'  ■"  "        '" 

,    L    (t)N    =    X 

I    y 
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Each    of    the    above    steps    requires    a    5-diagonal    matrix 

solver  which   will   be    described   in    the   Appendix  A. 

1. 2   On    the   boundary   points    of    the    computation    domain 

The    artificial    radiation   boundary    condition    is    of 
the    general    form 

*t  "^    ^'''x  +    V(t,j    =    0 
We    approximate    it  by 

^n+1    ,n-l  ^n+1^    n-1 

which    can  be    expressed    as 

(1+AtuD    +AtvD    )  d)N    =     (-1+AtuD    +AtvD    )  d)M 
X  y    ^  x  y    ^ 

-    2At(uD    +vD    ) *" 
X         y    ^ 

which  can  be  factored  within  first  order  in  space  and 

second  order  time  as 

(1+AtuD  )  (1+AtvD  )(!)N  =  RHS 
X        y  ^ 

The  algorithm  consists  of  the  following  two  steps. 
(1+AtuD  )X  =  RHS 

X 

(1+AtvD    )  ())N    =    X 

1 .  3   V7ake    condition 

As   before,   we    assume    that  both   pressure    and   normal 
velocity   components    are    continuous    across    the  wake  which 


59 


is    assumed    to    lie    on    a   segment   of    the    x-axis.       The  wake 
condition    is 

[(J)    ]    +    u[<t)    ]    =   0    on    the  wake 

where    u   is    the    average    velocity    of    the    upper    and    lower 
wake    velocities. 

The    finite    difference    approximation    is    given  by 

T     ,     ^         uAt       .1, 
Let   B   =   -rr-     then 
At 

t*i^    =  l.B  ''' 

Hence,    once    the    jump    at    the    tail   point   has    been    estimated, 
all    the    jump   in    the  wake    can  be    calculated   from   the   equa- 
tion   (1)  . 

We    remark    that    the    artificial    viscosity  we   have    added 
in    the    finite    difference    scheme    is    of    amount 

{h  (sign    u)  uu    ^   +   k  (sign    v)  vv      } 

+    y{h(sign    ")  u  (uu^^+vv^^)    +   k  (sign    v)  v  (uu^^+vv^^ )  } 

with    y    =    max    {0;     (1    -    —^)},h   =    Ax,    and  k    =    Ay.  '■  ' 

n 

The  term  in  first  braces  is  an  advection  viscosity  which  will 
damp  out  some  noise  generated  either  by  the  artificial 
boundaries  or  the  body  surface.   The  term  in  second  braces 
is  the  desired  shock  viscosity.   The  whole  artificial  vis- 
cosity can  be  cast  into  the   divergence  form 
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P       +    Q 

^  2       2 

with   P   =    h    {(sign    u)uu      +    p(sign  u) u (^    ^^ — )     } 

2_    2 
and   Q      =    k    {(sign    v)vv      +    y(sign   v)v("   ^      )     } 


Failure    to   maintain    proper    conservation    form   can    result 
in    computed    shock    speeds    that    depend   on    grid    spacing. 

2 .    Analysis    for  the    Finite    Differencing   Strategy    and 
The    Approximate    Factorization    Process 

VJe   have   shown    that    the    disturbance    information   in    the 
potential    flow    field   is    propagated    as    the    Doppler   sound  wave 
which    consists   of    the   advection    and  wave   propagation   effects, 
The   potential    flow   equation   is   nonlinear.      As    a   guide    to    the 
stability   of   the    difference    scheme,   we    consider   two    linear 
models,    the    advection    and    the  wave   equations. 

2.1    The    radiation   boundary    conditions    is    modeled  by    the    two- 
dimensional    advection   equation 

(j)^    +    u4)      +    vi})  .    =    0  (1) 

t  X  J 

Our    finite    differencing   strategy   says    that 

An+l_^n-l  ^        ^  n+1      n-1 

+     (uD    +vD    )  {^ ^ )    =0  (2) 


2fit  '      X        y'  '  2 

VJe   examine    the    amplification   of   a   Fourier   mode.       Substi- 
tuting   ())    =    (f)    e     ■         ■'    for    <t>    at    the    ic- level,    the    growing 
factor   (j)    is    governed  by 

*^    -    1    =   (pu(e'^^-l)  +   qv(e"^'^-l))  ((j)^+l) 
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for    the    case    that    u    >    0    and   v   >    0      where    p   =    At/Ax, 
q    =    At/ Ay,    E,    -   mAx,    n    =    nAy.       Hence, 


^ 


l-pu(l-cos    g)-qv(I-cos    ri)-i(pu   sin    g^+gv  sin    n) 
l+pu(l-cos    f;)+qv)l-cos    n)+i{pu   sin    C+qv  sin    n) 


<    1 


So   it   is    unconditionally    stable    for   this    case. 

The    scheme    for    the    other    cases   where    either    u,    or   v, 
or  both    of    then  may   be    negative    are    easily    shown    to  be 
unconditionally   stable. 

Next,    we    examine    the    approximate    factorization    method 
for    this    finite    difference    approximate    for    the    advection 
equation.       Our   approximate    factorization    says    that    (2)     can 
be    factored    as 


(1+    tuD^)  (1+    tvS    )(j)N    = 
X  y 


-  (1-AtuD    )  (1-AtvD    )  (})M    -2At(uD    +vD    \^ 
X  y    ^  X        y    ^ 


•>_.■._._       ,  ?K    i  (mx+ny) 

By    Fourier    analysis,    we    substitute    ())    =    4)    e 


/>,  *-i        r\ 


[1   +   pud-e""^^)  ]  [1   +   qv(l-e"^^)  ]  (({)^-4)) 
=   -[1    -    pu(l-e"^^)  ]  [1    -    qvd-e'^'^)  ]  {4)-l) 


-2  [pu(l-e"^^)    +    qvd-e"^"^)  ]({) 


* 


l-pu(l-cos    O+ipu   sin  g 
l+pu(l-cos    C)~ipu   sin  5 


l-qv(l-cos    n)+iqv  sin    x\ 
l+qv(l-cos    n)-iqv  sin    v\ 


<    1 


We,  therefore,  conclude  that  the  approximate  factorization 
method  does  preserve  the  unconditional  stability  of  our  finite 
differencing  strategy  and  that  our  numerical  scheme  for  the 
advection  equation  is  unconditionally  stable. 
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2.2   As  a  guide  to  the  difference  scheme  at  the  interior 

points  of  the  computational  domain,  we  consider  the  wave  equation 


tt     XX     yy 


Our  finite  differencing  yields 


,n+l  ,  _^n-l 
D..0>''  -  ^D_  4-  D,„)  *   -+$____ 


or 


tt^      XX    yy      2 


2  2 

[1  -  -^  (D   +D   )](j)N  =  [1  -  ^-(D   +D   )](|)M  +  At^lD   +D  )  ^^ 
2     XX   yy  ^  2    xx   yy   ^  xx   yy  ^ 


r^    ^    J.    J.    J.  ■  X  i(kt+mx+ny)    -,  n  .  .         vaj- 

SuJDStituting  ^   -    e  ^       and  letting  oj  =  i<At  , 

we  have 


2  2 

(cos  w-1)  =  [p  icos  5-  1)  +  q  (cos  n-  D]  cos  to 


cos  to 


2  2 

1  +  p  (1-  cos  C)+  q  (1-  cos  n) 


As  long  as  p,q,  are  real,  w  is  real  for  ail  ^  and  n  .   This 
means  that  the  finite  difference  approximation  is  unconditionally 
stable  . 
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Next,  our  approximate  factorization  preserves  this 
property  and  permits  us  to  solve  the  large  algebraic 
system  easily.   Indeed,  if  we  write 

2  2 

[  1  -  ^^  n  ]  [  1  -  ^^  D  ]  (i)N 
^     2   xx' '     2   yy 

2  2 

=  [  1  -  ^^  D   ]  [  1  -  ^^  D   ]  (})M  +  At"^  (D   +D   )  $" 
^      2    xx^ ^      2    yy  XX   yy  ^ 

-.  ^  ^     i  (kt+mx+ny)      , 
Let  4)  =  e  ,  we  have 

2  2 

1+p  q  (1-cos  E,)   ( 1-cos  n) 

cos  oj  -   y * — 2 2 — 2  ~ 

1+p  (1-cos  ^)+q  (1-cos  n)+p  q  (1-cos  C)  (1-cos  n) 

For  all  real  p  and  q,  lo  is  real  if  C  and  n  are  real.   In 

other  words,  t  is  real  whenever  x  and  y  are  real.   This  means 

that  the  scheme  is  unconditionally  stable. 

Finally,  we  remark  that  the  scheme  has  no  time  step 

At  restriction  based  on  a  linear  stability  analysis. 

However,  in  actual  computation,  an  instability  can  be 

generated  by  the  motion  of  shocks  across  which  the 

differencing  switches  from  upwind  to  central.   To  prevent 

this  instability  from  occurring,  it  has  been  found  in 

practice  that  the  time  step  At  must  be  chosen  small  enough 

that  such  shocks  do  not  move  a  distance  greater  than  one 

spatial  grid  point  per  time  step.   This  restriction  is 

necessary  to  maintain  time  accuracy  anyway,  and  it  is  much 

less  severe  than  the  time  step  At  restrictions  associated 

with  explicit  methods. 
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V.   COMPUTATIONAL  RESULTS 


Our  computer  code  UFL05  consists  of  steady  and  unsteady 
modes.   The  steady  mode  is  the  standard  line  relaxation 
scheme   for  the   steady  equation.   We  use   it  to  generate 
a  good  initial   guess   for  the   unsteady  mode.   In  fact, 
any   steady  potential  flow   solver   can  be   used  to 
replace  this  steady  routine.   The  unsteady  mode   can 
also  be  used  to  compute  the  steady  solution.   In  this  section 
we  first  check  the  unsteady  mode  by  calculating  some  steady 
solutions.   Then  we  present  some  computational  results  for 
conventional  and  supercritical  wing  sections  in  rigid  body 
motion. 

1.  Steady  Calculations 

As  a  test  case,  steady  state  calculations  for  the 
NACA0012  airfoil  at  Mach  number  M  =  0.79  and  angle  of 

oo 

attack    a  =   0°    are   performed  by    the   standard    line    relaxa- 
tion  method    for    the    steady    equation    and  by    the    unsteady 
scheme.      The    two   modes    produce    virtually   identical    results. 
The    time    step   size    for   the    unsteady   mode    in    this    calculation 
is    set    to    At   =    lOAx  which    is    much    larger    than    the    time    step 
allowed  by    the    CFL    condition    for    the    explicit   method.       For 
a   coarse   mesh   of    32    x    8    grid  points,    the    time    required   to 
converge    to    the   steady   state    using   the    unsteady   mode    is 
comparable    to   the   steady   mode. 
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The  optimal  location  of  the  artificial  boundary  is 
problem  dependent.   If  the  artificial  boundary  is  moved 
too  close  to  the  airfoil,  instability  can  occur.   The 
computational  domain  shown  in  Figure   8c  has  also  been 
used  for  the  NACA0012  airfoil  at  M  =0.79  and  a  =  0°. 

00 

Note  that  the  upstream  boundary  is  about  1.5  chord 
lengths  from  the  nose  as  compared  to  a  distance  of  about 
10.5  chord  lengths  used  for  the  above  example.   The  ratio 
At/Ax  is  given  the  value  10  as  above.   The  correction  is 
observed  to  decrease  much  more  slowy  in  this  case.  Since  the 
grid  system  is  stretched  in  the  code,  the  reduction  in  the 
computational  mesh  is  not  linearly  proportional  to  the  physical 
distance  of  the  boundary  from  the  airfoil.   The  benefit  obtained 
by  the  reduced  number  of  mesh  points  is  overshadowed  by  the 
reduced  numerical  stability. 

There  is  no  difficulty  in  calculating  flows  with  sonic 
flight  speed.   A   Joukowski   airfoil  at  M^  =  1  and  a  =  0° 
is  chosen  as  an  example,  v;ith  the  ratio  At/ Ax  set  to  3.5 
in  this  case.   Usually,  the  numerical  stability  of  the  un- 
steady mode,  in  terms  of  the  ratio  At/Ax,  decreases  with 
either  flight  speed  or  angle  of  attack.       _^  . 

2 .  Unsteady  Calculations  . 

The  NACA0012  airfoil  and  KORN  airfoil  (75-06-12)  are 
chosen  as  prototypes  for  conventional  and  supercritical 
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airfoils,  respectively.   Botli  airfoils  have  the  same  thick- 
ness to  chord  ratio.   The  rigid  body  motion  of  an  airfoil 
can  be  described  by  three  parameters:  angle  of  attack, 
flight  speed,  and  flight  angle.   We  consider  the  flow  past 
each  airfoil  when  these  parameters  are  varied  separately. 
2 . 1  Variation  of  flight  speed 

First,  we  consider  the  acceleration  of  the  airfoil  in 
the  streamwise  direction.   That  is,  the  airfoil  moves  with 
a  sinusoidal  variation  in  flight  speed  but  with  flight 
angle  and  angle  of  attack  fixed. 

In  Figure  11  a,D,c  ,  for  the  NACA0012  airfoil,  as  the  flight 
speed  increases  (decreases),  the  supersonic  region  grows 
(shrinks)  in  size  and  the  shock  strengthens  (weakens)  and 
moves  aft  (fore).   The  shock  wave  displacement  can  also  be 
observed  in  the  pressure  distributions  in  Figure  11  d.e.f.,  or 
from  the  position  of  peaks  in  the  traces  of  pressure  sensors 
on  the  upper  surface  of  the  airfoil  in  Figure  11  g  .   The 
peaks  in  those  local  pressure  traces  are  produced  as  the 
shock  wave  passing  by  the  pressure  sensors.   The  nonsinusoidal 
trace  curves  demonstrate  the  nonlinear  transonic  effects 
caused  by  the  shock  wave  displacement.  The  same  calculations 
for  the  KORN  airfoil  appear  in  Figure  12.   The  unsteady 
loading  distributions  are  shown  in  Figure  12  e,f,     where 
peaks  in  the  loading  distributions  are  again  due  to  shock  waves, 
The  loading  is  the  difference  of  lower  and  upper  pressure 
coefficients . 
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2 . 2  Variation  of  angle  of  attack 

Next  we  consider  the  pitching  airfoil  which  moves  with 
a  sinusoidal  variation  of  the  angle  of  attack,  but  with  the 
flight  velocity  fixed.   Small  variations  in  angle  of  attack 
may  lead  to  considerable  changes  in  the  pressure  distribu- 
tion, shock  position  and  shock  strength.   As  shown  in 
Figure  9  a,b/r/  for  the  NACA0012  airfoil,  when  the  angle  of 
attack  increases  (decreases),  the  supersonic  region  on  the 
upper  surface  of  the  airfoil  grows  (shrinks)  in  size,  and 
the  shock  wave  strengthens  (weakens)  and  moves  aft  (fore). 
In  Figure  9j  the  nonsinusoidal  trace  of  the  pressure  at  loca- 
tion *6  on  the  airfoil  surface  clearly  displays  the  shock  wave 
movement.   The  unsteady  pressure  loading  distributions  are 
shown  in  Figure  9g,h,i. 

The  same  calculations  were  performed  for  the  KORN  airfoil 
as  shown  in  Figure  10.    It  is  worthwhile  noticing  that  the 
unsteady  traces  of  the  local  pre-sure  sensors  for  the  KORN 
airfoil  are  much  more  nonlinear  than  those  for  the  NACA0012 
airfoil.  This  pattern  is  also  observed  in  the  loading  distri- 
bution for  the  two  airfoils.   The  fact  that  the  shock 
excursion  amplitude  decreases  with  an  increase  in  oscillatory 
frequency  can  be  seen  from  the  unsteady  traces  of  the  pressure 
sensor  *6  in  Figure   9  j,k,l. 

2.3  Variation  of  flight  angle 

Finally,  we  consider  changes  in  the  airfoil's  flight 
angle  while  keeping  ±he   angle  of  attack  and  flight  speed 
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fixed.       The   motion,    for   angle    of   attack   a  =   0,    of    the    air- 
foil  is    described  in   Figure    7. 


Figure    7.       The   Chords   of   Airfoils    are    Tangent    to    the   Flight 
Trajectory. 

The   characteristics    of    this   case    are   similar   to    those   of 

the   pitching   airfoil.      Computational    results    for    the 

NACA0012    and  K0R1^I    airfoils    are    shown   in   Figures    13      and 

14     ,    respectively. 
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VI.   CONCLUSION 

1 .  Summary  of  the  Work 

A  numerical  method  has  been  presented  for  determining 
the  inviscid  transonic  flow  past  airfoils  in  rigid  body 
motion.   The  method  is  based  on  the  unsteady  transonic  po- 
tential flow  equation  in  a  computational  domain  designed 
for  accurate  application  of  the  body  surface  boundary  con- 
dition. A  set  of  first  order  radiation  boundary  conditions  are 
applied  at  the  artificial  computational  boundaries  located  at  a 
finite  distance  away  from  the  airfoil  surface.  The  finite  differ 
ence  approximations  of  the  potential  flow  equation  and  radiation 
boundary  conditions  are  constructed  by  using  a  type  dependent 
differencing  strategy.  The  leading  truncation  term  provides  the 
necessary  dissipation  to  stabilize  the  scheme  and  to  capture  the 
shock  waves  automatically.   The  finite  difference  approximations 
are  perturbed  within  the  same  order  of  accuracy  in  order  to 
permit  their  factorzation  into  one-dimensional  operators.  Conse- 
quently, the  problem  can  be  solved  efficiently  by  using  a 
5-diagonal  matrix  solver.  The  resulting  algorithm  is  a  time 
marching  scheme  without  any  iteration  process  in  each  time  step. 

Numerical  experiments  show  that  the  scheme  is  very 
stable.   The  results  presented  in  Section  V  demonstrate 
that  the  method  is  able  to  resolve  the  highly  nonlinear 
transonic  flow  effects  for  flutter  analysis  of  airfoils  as 
long  as  the  boundary  layer  remains  attached. 
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2 .    Extension    of    the    Technique 

Our  numerical   method   can  be   extended   to   three    space 
dimensional   problems.      An    important    application    is    the    un- 
steady   transonic    flow   past  wing-body    combinations    that 
model   our   airplane.       The   necessary    geometric   mapping   tech- 
niques   are    available    in    analysis    codes    that   compute    steady 
flow   past    a  wing-body    combination    [27].       Singularities 
associated  with    the    geometric   mapping  would   not  be    a 
serious    problem   and   could  be    treated    in    a   manner   similar 
to    the    one    used   in    the    steady    calculation.      A    further 
application   could  be    the   helicopter   rotor   in    forward 
flight.      Here    the    flow   is    unsteady   because   of    the    relative 
velocity   of    the    advancing    and   retreating  blades     [8], 

In   principle,    shock    accuracy  would  be    improved  by 
shock    fitting  methods    [45]    or,    alternatively,    by    the    use 
of   a   difference   scheme   in    conservative    form.       In   practice, 
however,    a    turbulent  boundary    layer    correction   would  be 
needed   for   more   exact  shock    jump   modeling    [10,    13,    35], 
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APPENDIX 

A .    A  5-Diagonal  Matrix  Solver 

Here,  we  present  a  method  of  solving  a   5-diagonal 
matrix  problem.   Suppose  Ax  =  y  is  solved  for  x,  where 
X  and  y  are  nxl  column  vectors  and  A  is  of  the  form 


^1   ^1   ^1 

^2   ^2   ^2   ^2 

d3   b3   33   C3   e3 


n-2  n-2  n-z  n-z  n-2 

d  -,  h     ,  a  ,  c  ^ 
n-1  n-1  n-1  n-1 

d   b   a 

n    n   n 


Assume  the  matrix  can  be  factored  in  the  tridiagonal  form 


A  =  LU  = 


ct-, 


a. 


n 


n   n 


1   Yi   e^ 


^2   "2 


1  Y  -  e  - 
'n-2  n-2 


1  Y, 


n-1 
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Then  we  find 

6.  =  d. 
1     1 

(1)         6.e._2  +  B.Yi_i  +  "i   =   ^i 

B.e.  ,  +a.Y-   =   c. 
1  1-1     I'l      1 

a  .  =  e . 
1    1 

for  i  =  l,...,n   if  the  default  values  are  set  equal  to 

zero.   Namely,  b,  =d,  =d„=c   =e   ,=e   =0, 
-'I     1     2     n    n-1     n 

6,  =6,=6„=Y   =e   i=e   =0   and  which  can  be 
1     i     2     n    n-i     n 

solved  as 


6.  =  d. 
1     1 


(2)  a.  =  a.  -  B-Y-  i  -  5.e.  ^ 

1     1     1 ' 1-1     1  1-2 


Yi  =  (c.-  6ie._^)/a. 


e  .  =  e./ot. 
1     11 


for  i  =  l,...,n,  in  ascending  order  if  none  of  the  a.  vanish. 
The  intermediate  step   Lg  =  y  becomes 

6^g^_2  +  3j^g^_-L  +  a^^g^^  =  Y^   for   i  =  l,...,n, 

we  can  solve  this  system  recursively  in  ascending  order. 
Namely, 


(3)       g^  =  (yi-Bigi.i-5igi_2)/ai 


for  i  =  l,...,n,  if  the  default  values  g_-,  =  g_2  =  0- 
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The  final  step  Ux  =  g   is   expressed  by 


X. 


i  -"   ^i^i+1  ■"  ^i^i  +  2  =  ^i   ^°^   i  =  l,...,n. 


where   x  ,,  =  x  ,„  =  0.   The  system  can  be  solved  in 
n+i    n+Z  ■" 

descending  order   recursively  as 


X.  =  g.  -  YiX.^i  -  eiX.^2   ^°^   i  =  n,...,l. 


We  remark  that  the  LU  factored  form  is  not  unique; 
for  example,  L  and  U  can  be  the  following 


L  = 


1 

63   B2   1 


n   n 


U  = 


°'l   "^1   ^1 


"2   ^2   ^2 


a  ~  Y  ->  E  o 
n-2   'n-2   n-2 


ct   1  Y   n 
n-1  '  n-1 


a 


n 
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B .    Computer  Program  UFL05 
1.    Operation  of  the  Program 

The  sheared  parabolic  coordinates  described  in  Sec- 
tion 111,3   are  introduced.   The  input  parameters  XSING 
and  YSING  determine  the  location  of  the  singular  point 
about  which  the  square  root  transformation  is  made.  It 
is  important  to  choose  these  two  parameters  so  that  the 
unfolded  profile  does  not  have  any  sharp  bumps.   The 
mapped  coordinates  are  printed  so  that  this  can  be 
checked. 

The  difference  scheme  for  the  steady  routine  used  to 
initialize  the  calculation   is  in  fact  the  standard  line 
relaxation  method.   Faster  convergence  is  usually  obtained 
by  using  horizontal  relaxation,  y-sweep,  marching  toward 
the  body.   The  difference  scheme  for  the  unsteady  routine 
conforms  closely  to  the  description  in  Section  IV,  1.   It 
is  implemented  in  the  computational  domain  described  in 
Section  III,  3   as  first  performing  a  y-sweep,  marching 
toward  the  body  with  horizontal  lines,  then  followed  by  an 
x-sweep,  with  vertical  lines  marching  from  left  boundary 
toward  right  boundary   of  the  computational  grid. 

The  initial  values  of  the  time  dependent  problem  are 
provided  by  either  using  unsteady  mode  alone  or  using  both 
steady  and  unsteady  modes.   The  program  contains  a  switch 
for  the  choice.   For  fine  mesh,   such   as   128  x  32,  the 
method  employing  both  modes  is  recommended.   A  run  using 
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both  modes   can  be  described  as  follows.  First,  using 
steady  mode,  calculations   are  first  performed  on  a 
coarse  mesh   and  then  on  a  fine  mesh  with  twice  as 
many  intervals  in  each  coordinate  direction.   The  coarse 
mesh  result   is   interpolated   to   provide  the   starting 
guess   for  the  fine  mesh.   It  usually  consists  of  200 
cycles  on  coarse  mesh,  32  x  8,  followed  by  100  cycles  on 
a  fine  mesh,  64x16,   50  cyles  on  a  finer  mesh   128^32. 
The  resulting  reduced  velocity  potential  is  used  as  the 
starting  guess  for  the  steady  iteration  using  the  unsteady 
routine.   After  75  cycles  in  this  steady  iteration  step, 
we  begin  our  time  marching  calculation.   A  better  initial 
value  can  be  obtained  after  one  complete  periodic  cycle. 
Computational  results  show  that  the  difference  between  the 
second  and  the  third  period  cycles  is  small.   We  therefore 
consider  the  results  from  the  second  period  as  our  desired 
output. 

The  input  data  deck  for  a  run  is  arranged  to  include 
title  cards  listing  the  required  data  items.   The  complete 
set  of  title  cards  provides  a  list  of  all  the  data  which 
must  be  supplied  and  can  be  used  as  a  guide  in  setting  up 
the  data  deck.   Each  title  card  is  followed  by  a  card 
supplying  the  numerical  values  for  the  parameters  listed. 
The  input  parameters  are  given  in  the  Glossary  in  the  order 
of  their  appearance  on  the  data  cards.   All  data  items  are 
read  in  as  floating  point  numbers   in  fields  of  10  columns. 


and  values  representing  integer  parameters  are  converted 
inside  the  program.   The  data  deck  for  NACA  0012  at 
M=0.79,   a=0°  +  l°  sin  kt  is  shown  in  Table  1  . 
The  output   consists   of  printout   and  Calcomp   plots. 
The  program  prints  the  mapped  coordinates  of  the  airfoil 
generated  at  the  mesh  points  of  the  computational  grid. 
Parameters  such  as  mesh  size,  flight  speed,  flight  angle, 
angle  of  attack  are   also  printed  so  that  the  case  can  be 
identified  easily. 

For  each   iteration  using  the  steady  routine  the 
program  prints  the  iteration  number,  the  maximum  correction 
to  the  reduced  velocity  potential,  and  the  maximum  residual 
for  the  steady  flow  equation   together  with  the  coordinates 
of  the  point  where  these  occur  in  the  computational  grid,  the 
circulation,  the  relaxation  factor   pi,  p2 ,  p3 ,  and  the 
number  of  supersonic  points.   After  a  maximum  number  of 
cycles  has  been  completed  or  a  convergence  criterion  has  been 
satisfied,  the  angle  of  attack,  flight  speed,  flight  angle, 
lift,  drag   and  moment  coefficients   are  printed.   If 
desired,  the  pressure  distribution  along  the  airfoil  surface 
and  a  chart  of  the  local  Mach   numbers  can  be  printed. 
If  the  mesh   is   to   be   refined,  the  program  then  repeats 
the  same  sequences  of  calculations  and  output  on  the  same 
mesh.   A  Calcomp  plot  is  generated  to  show  the  pressure 
distribution  over  the  airfoil  on  the  finest  mesh  at  the 
end  of  this  subroutine. 
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For  a  steady  iteration  using  nonsteady  mode,  the 
program   first  prints  the  flight  conditions,  the  mesh  size, 
and  the  dimensionless  time  step.   After  each  iteration, 
the  program  prints  the  maximum  change  in  the  velocity 
potential  with  the  coordinates  of  the  point  in  the  grid 
system.   If  desired,  the  pressure  distribution  along  the 
airfoil  and  the  local  Mach  number  chart  can  be  printed. 
Calcomp  plots  for  the  pressure  distribution,  the  leading 
distribution,  and  the  supersonic  zone  over  the  airfoil  are 
generated  separately  at  the  end  of  this  step. 

Before  the  unsteady  time  marching  process,  the 
advanced  time  steps  required  to  finish  the  assigned  period 
is  estimated  and  printed.   After  one  complete  period  has 
been  computed  the  flight  conditions   together  with  the 
aerodynamic  forces,  lift,  drag  and  moment  coefficients,  are 
thereafter  printed  periodically.   If  desired,  the  pressure 
distribution  over  the  airfoil  is  also  printed.   Calcomp 
plots  for  the  pressure  distribution,  the  loading  distribu- 
tion,  and  the  supersonic   zone  over   the  airfoil   are 
generated  periodically.   At  the  end  of  the  calculation  the 
unsteady  traces  of  the  airfoil  motion  and  the  grid  system 
near  the  airfoil  are  plotted.   The  graphs  can  also  be 
produced  as  individual  frames  in  a  film  strip.   Then  a 
complete  history  of  the  time  dependent  motion  will  be  visible. 
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2.      Glossary  of  the  Program 

The  input  parameters  are  lifted  in  the  order  of 
their  occurrence  on  the  data  title  cards. 

Title 
Card  1 

ISYM    Indicates  the  type  of  profile. 

ISYM  =  0  denotes  a  cambered  profile.   Coordinates 

are  supplied  for  upper  and  lower  surfaces,  each 

ordered  from  nose  to  tail   with  the  leading  edge 

included  in  both  surfaces. 

ISYM  =  1  denotes  a  symmetric  profile. 

A  table  of  coordinates  is  read  for  the  upper  surface 

only. 
NU     The  number  of  upper  surface  coordinates. 
NL     The  number  of  lower  surface  coordinates. 

For  ISYM  =  1,  NL  =  NU   even  though  no  lower  surface 

coordinates  are  given. 
NX     The  number  of  mesh  cells  in  the  direction  of  the 

chord  used  at  the  start  of  the  calculation.  NX  =  0 

causes  termination  of  the  program. 
Ny     The  number  of  mesh  cells  in  the  direction  normal  to 

the  chord. 
MHALF   Determines  whether  the  mesh  will  be  refined. 

MHALF  =  0.  The   computation  terminates  after  completing 

the  prescribed   number  of  iteration  cycles  or  after 

convergence  for  the  input  mesh  size. 
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MHALF  7^  0.   The  mesh  spacing  will  be  halved  after 
NRELAX  cycles  have  been  run  on  the  crude  mesh  size. 
The  refinement  will  be  performed  MHALF  times. 

RSTAD   Determines  whether  the  steady  mode  will  be  employed. 
RSTAD  =  0.   The  steady  mode  will  not  be  called, 
the  steady  flow  calculation  entirely  depends  upon 
the  unsteady  mode.   RSTAD  =  1.   Both  steady  and 
unsteady  modes  are  employed  for  the  steady  flow 
calculation. 

STADI   Indicates  the  type  of  flow  calculation. 

STADI  =  1.   The  computation  is  running  for  the 

steady  state  solution.   STADI  =  0.   The  computation 

is  a  time  dependent  run. 

Title 
Card  2 

NRELAX  The  maximum  number  of  iteration  cycles  which  will  be 
computed  in  the  steady  iteration  process. 

RELAXTO  The  desired  accuracy.   If  the  maximum  correction 
is  less  than  RELAX  TO   the  calculation  terminates 
or  proceeds  to  a  finer  mesh;  otherwise  the  number  of 
cycles  set  by  NRELAX  are  completed. 

CHECKPT  Determines  whether  the  CHEKPTX   is  required. 

CHEKPTX  =  1.   The  CHEKPTX  is  called. 

Title 
Card  3 

COORS   The  stretching  factor  in  the  x  coordinate  stretching 

transformation  described  in  Section  III,  3. 
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COORT  The  stretching  factor  in  the  y  coordinate  stretching 

mapping  described  in  Section  III,  3. 

RCBDY   To  locate  the  computational  boundaries. 

Title 
Card  4 

PIO     The  subsonic  relaxation  factor  for  the  reduced 

potential  in  the  steady  flow  calculation  routines. 

It  is  between  1.  and  2.   and  should  be  increased 

towards  2.  as  the  mesh  is  refined. 
P20     The  supersonic  relaxation  factor  for  the  reduced 

potential  in  the  steady  routine.   It  is  not  greater 

than  1.  and  normally  set  to  1. 
P30     The  relaxation  factor  for  the  circulation.   It  is 

usually  set  to  1.   but  can  be  increased 

PlOl 

P102    The  increments  of  PlO  as  the  mesh  system  is  refined 

P103 


1  time,  2  times  and  3  times,  respectively, 


Title 
Card  5 


FRFORA 

j^j^pi^-   The  frequency  rate  (rad/time)  and  amplitude  (degree) 

of  the  sinusoidal  variation  of  angle  of  attack  (in 

degrees) . 

FREORM 

AMPLM  The  frequency  rate  and  amplitude  (mach  number)  of  the 

sinusoidal  variation  of  flight  speed  in  mach  number. 

FREQRC 

AMPLC   The  frequency  rate  and  amplitude  of  the  sinusoidal 

variation  of  flight  angle  in  degrees. 
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PERIOD  The  complete  sinusoidal  periods  to  be  calculated. 
DEGREE  The  degree  interval  to  plot  the  graphs,  pressure 

distribution,  the  loading  and  the  supersonic  zone 

over  the  airfoil. 

Title 
Card  6 

ALPHAl  The  angle  of  attack  in  degrees. 

MACHl   The  flight  speed  in  mach  number.  The  speed  of  sound 

at  infinity  is  set  to  be  unity. 

THETAl  The  flight  angle  in  degrees 

TSRATIO     The  ratio  At/Ax,  .  v   between  time  step  and 

(min)  '^ 


minimum  spacial  step. 


Title 
Card  7 


TE  ANGLE    The  included  angle  at  the  trailing  edge  in  degrees. 
The  profile  may  be  open,  in  which  case  it  is  the 
difference  in  angle  between  the  upper  and  lower 
surfaces. 

TE  SLOPE    The  slope  of  the  mean  camber  line  at  the  trailing 

edge.   This  is  used  to  continue  the  coordinate  surface, 
assumed  to  contain  the  vortex  sheet,  smoothly  off  the 
trailing  edge. 

Y  C  T  \T^^ 

YSTNC  '^^^   coordinates  of  the  singular  point  inside  the  nose 
about  which  the  square  root  transformation  is  applied 
to  generate  parabolic  coordinates.   This  point  should 
be  located  as  symmetrically  as  possible  between  the 
upper  and  lower  surfaces  at  a  distance  from  the  nose 


143 


roughly  proportional  to  the  leading  edge  radius. 

It  can  be  seen  whether  the  location  has  been  correctly 

chosen  by  inspecting  the  coordinates  of  the  mapped 

profile  printed  in  the  output.   If  the  mapped 

profile  has  a  bump  at  the  center,  the  singular  point 

should  be  moved  closer  to  the  leading  edge.   If  the 

mapped  profile  is  not  symmetric  near  the  center, 

with  a  step  increase  in  y,  say,  as  x  increases  through 

0,  the  singular  point  should  be  moved  closer  to  the 

upper  surface. 

Title 
Card  8 

X 

Y  The  coordinates,  upper  surface  coordinates,  of  the 
THETA 

upper   surface   and  its  tangent  angle  in  degrees. 

These  are  read  on  the  data  cards  which  follow,  one 
pair  of  coordinates  and  its  tangent  angle  per  card 
in  the  first  three  fields  of  10  from  leading  to  trail- 
ing edge  inclusive. 

Title  '  _ 

Card  9 

X 

Y  The  coordinates  and   its  tangent  angle   at  the  lower 
THETA 

surface,  read  from  leading  to  trailing  edge.  The 

leading  edge  point  is  the  same  as  the  upper  surface 
leading  edge  point.   The  trailing  edge  point  may  be 

different  if  the  profile  has  an  open  tail. 

Title 

Card  10  ....... 

End  of  the  calculation. 
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"''"^'s,^^    Cols. 
Card       ^^-^v,^^ 

1-10 

11-20 

21-3C 

31-40 

41-50 

51-60 

61-70 

71-80 

NACA0012 

1 

ISYM 

NU 

NL 

NX 

NY 

MHALF 

RSTAD 

STADI 

1. 

37. 

37. 

32. 

8. 

2. 

1. 

0. 

2 

NRELA) 

RELAX 
TO 

CHECKPT 

' 

200. 

l.E-6 

0. 

3 

COORS 

COORT 

RCBDY 

0. 

1. 

1. 

4 

PIO 

P20 

P30         PlOl 

P102 

P103 

0.94 

0.8 

1.               .19 

1 

.58 

.72 

5 

FRBQRA 

AMPLA 

FREO^ 

AMPijyi 

FREQC 

AMPTC 

PERIOD 

DEGREE 

0.3            1. 

0. 

0. 

0. 

0. 

2. 

90. 

6 

ALPHAl 

mach: 

IHETAl 

TS RATIO 

' 

0. 

.79 

0. 

10. 

7 

TEANGLE   I'ESLOPE 

XSING 

YSING 

16.15         0. 

.8 

0. 

8 

X                 Y 

THETA 

(UPPER   SURFACE) 

1 
NN 

■ 

9 

X 

Y 

THETA 

(LOWER   SURFACE) 

1 
NL 

Table  1.   Data  Deck  for  the  Program. 


145 


LISTING  GF  THC  PkUGRAM 


PROGRAM 

i 
THE  ANA 
THE  UnS 
LiGUNUAi^ 
ARE  SQL 
FIVE  01 
PPGGPAM 
G  IS  TH 
CuMMUN/ 

1 

2 

3 
CCMMUN/ 
CCMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 

1  ,AMPLM 

?  ,FRFwR 
COMMON/ 
COMMON/ 

1 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 
DATA  VA 
IkEAL;  = 
I  U  R  I  T 
IPLOT 
PI 
PAD 
.  WRITE  ( 
V.  R  I  T  E  ( 
I     FORMAK 

1       5 


0?LDt> 

TAPE 

LYSIS 

TEAUY 

Y  CGND 

VEU  BY 

AGONAL 

MtD  BY 

E  VELC 

A/  GM( 

,  A0( 

,63{ 

*  AL, 

3/  SVC 


C/ 

D/ 
E/ 
F/ 
&/ 
H/ 
1/ 
J/ 


XP( 

SLO 
CHO 
XR> 
TIT 
DX, 
X(2 
RAU 


(  INPOT 
7) 

;;f  tra 

TPANSU 

ITIONS 

AN  AL 

MATRI 

I-CHU 

CITY  P 

132>  3o 

132)  »A 

36) ^NX 

UTIM,C 

132)»S 

26C) >Y 

PT>TRA 

RD*XM, 

yr,ks^ 

L  E  (  2  0  ) 

Dt,DT, 
feO)> Y( 
*PI^  AL 
M,FASA 
E 

Il>12* 
(801), 
*  N  I  T  S  , 
P2,P3, 
RS^COO 
■?R»  IP 
NIT, WG 


, OOF  POT, TAPE  5  =  INPUT* TAP Eb  =  OUT PUT* TAPE98«0UTPUT, 


NSQNIC 
NIC  POT 
IN  MOV 
TERN  AT  I 
X  SOLVE 
NG  CHAN 
OTENTI A 
),G( 132 
1  (132), 
,NY,  IXl 
B,Sd,NS 
M  {  1 3  2  )  , 
P{260), 
IL,SCAL 
CL,CD,C 
XS(50C) 
,  IPLOT 
DXX, DYY 
260) 
S, ALT, A 
GM,CETA 


FLOW  PAST  AIRFOIL  IN  RIGI 
ENTIAL  FLOW  EQUATION  WITH 
ING  SHEARED  PARABOLIC  COO 
NG  DIRECTION  IMPLICIT  SCH 
R 

G  DURING  SEPTEMBER  198C 
L  IN  THE  ABSOLUTE  FRAME 
,36),GN(132,36), S0(132  ),S 
A2( 132),A3(132), 80(36)  ,B1 
, 1X2, KSYM,FM AC H, ALPHA, CA, 
,  R  G  ,  I  6  ,  J  G 
CP(132) 
Dl(260),D2(260),D3(2bO) 

M 
,YS(bOO) 

,DTT,DXY,OXT,DYT,TSR 

LTT,AMPLA,FREQRA,FASAGA,F 
S,CFTAT,CETATT,AMPLC,FREQ 


C  BODY  MOTION 

RADIATION 
RDINATE  SYSTEM 
EME  WITH 


1(132), S2(132) 

(36),B2(36) 

SA,FMACH2 


I3,J1,J2,J3 
TC0(B01),TCM(d01),TCP(9>801),TCPS(9),CLS,CUS 

1JUMP,NSTEP,JSTEP,PERI0D,MHALF 

TAU 

RT 

,  JR,  IRSTAD       '  ■  ' 

(132  )  .     '  ,    ■ 


READ 
WRITE 


, FPEQk 
,  IPSUP 
K/  10, 
L/  TCL 
,  C  M  S 
M/  PI, 
0/  COG 
STADI/ 
WAKE/ 
R/0/ 
5 

=  6 

=  -1 

=  3  .1^15926  53  589  79 

=  57.2957795130823 
IWR1T,600) 
IWRIT,2) 

l^HOPROGkAM  UF L J  5, 70 X,  32H  I-ChUNG 
6H0S0LUTI0N  OF  UNSTEADY  TRANSONIC 
IREAD,530)  ( TITLE( I ),I«1,20) 
IWRIT,630)  (TITLE (  I),I  =  1,20) 


MACHS,FMACHT 
RCFASAGCCETA 


CHAi>IG,COURANT    INSTITUTE/ 
POTENTIAL    FLOW    EQUATIONS) 
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PEAD 

READ 

ISYM 

IKSTAO 

IS  TAD  I 

NU 

NL 

IF  (NU 

NXO 

NYO 

NX 

NY 

IF(NX. 

MHALF 

NHALF 

READ 

REAL) 

ICHECK 

('E^D^  I 

Rt A0(  I 

LhALF= 

1F(LHA 

READ( I 

REALK  I 

KEAD(  I 

RtAD(  I 

IF (PER 

IF (PER 

DEGPE= 

I  F  (  P  E  R 

DEGRE= 

IGRAF* 

IF( IGR 

IPSURE 

FkEQR= 

IF (FkE 

FkEOR  = 

IF{FRE 

freqr» 

I F  (  F  R  E 

FREQP= 

IF(FkE 

MITO 

PEAD 

READ 

FMACh 

FMAChZ 

CETA  = 

CETAS= 

RtAD 

READ 

TRAIL 

N 

READ 


FSYM,FNU»FNL»FNXfFNY,FHALF,FRSTAD*FSTADI 


TO    3C2 


(  IPEAD,  tjOC) 
(  IREAD*  51C) 

=  F  SYM 
=  FRSTAD 
=  FSTADI 

=  FNU 

=  FNL 
.LT.l)  GO 

=  FNX 

=  FNY 

=  NXO 

=  NYO 
NE.^*NY)  GU 

=  FHALF 

=  0 
(  IRE  AD»  500  ) 
( IREAD, 510) 
=  CHECPT 
READ>50C) 
READ>510) 

RCBDY 
LF.LE.O .Ok, 
READ>50C) 

READi510)  P10fP20fP30tP101^P102»P103 
READ»500) 

PEAD»51C)FREQRA,AMPLA,FRE0kM,ANPLM,FR£QRC»AMPLC>PERIQD>0EG«E 
I0D.GT.2.)  DEGRt^  <.  5  . 

LT.3)  G'J    TO  3 


TO  3J2 

FIT,COV>CHECPT 

CDCRS>COQRT>RCBDY 
LMALI-.GT.3)  LHALF  =  1 


12 


IQD 

45. 
I0D.LT.4 )     GC    TO    3 

90. 

3(30./DEGk£ 
AF.GT.12)     IGRAr= 
=     PERICD+IGRAF 

ICO. 
QRA.LE.O.)    GU    TO    4 

Ar^INK  FRFQRA,FREOR) 
ORM.LE.O.)  GO  TO  5 

AMINK  FREQRM>  FREQR) 
QRC.LE.O.)  GO  TO  6 

AKINK  FREURC*  FREQR) 
QR.LE.O. )  GO  TO  302 

=  F  IT 
(IREAD,50C) 

(IREAD>5iO) 
=  FMl 
=FMACH*FMACH 
CTl 

CETA/RAD 
(  IREAD*  500 ) 

(IkEAD,510)  TRAiL*SLOPT/XR,YR 
=  TkAIL/RAD 
=  NL   +NU   -1 
(  IREAD*500 ) 


ALU  F"11,CT1,TSR 
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10 
11 


12 


13 


37 

36 


1^ 


16 


DU    7 

f-  EAU 

L 

IF     (I 

ktAD 

DlI    b 

READ 

J 

X(  J  ) 

Y(  J) 

GG    Tu 

J 

OL     10 

J 

X(  J) 

Y(  J  ) 

CHORD 

XM 

AL 

ALPHA 

ALS  = 

KSYM 

IF     {  A 

CALL 

IF(  IX 

I  r  (  I  K 

CALL 

bCj    TLj 

CALL 

UTIM  = 

MIT 

ALT  =  0 

ALTT  = 

FMACh 

FMACH 

Cf  TAT 

CETAT 

WFITE 

w  k  I T  e 

KX=  N 
DO  lb 
WRITE 
WRITE 
WRITE 
MY=  N 
DU  lb 
«  k  I  r  E 
IF  (  IR 
IF  (NH 
NIT=U 
Jl=  3 
KhALF 
J3=  3 
J2»    J 


-X(IIL) 

+.2b*CHURD 


I  =  N  L  ^  N 
( Ikb AD, 510 )     X(  1)>  Y(  I  ) 
=    NL       +1 
SYM.GT.O)     GO    TJ    9 

(  IRE  AD,  5)00) 
1=     IWiL 
(IkEADfblC)     VAL^UUfi 
=    L       -I 
=    VAL 
=     PUM 
11 

=    L 
I  =  ->!  L  ,  N 

=  J,   -1 
=  X(I  ) 
=  -Yd  ) 
=  X  (1) 
=  X  (NL) 
=  ALl 
=  AL/kAD 
ALPHA 

=  iSYM 
LPHA.NE.O. )  KSYM  =  0 
CQORDCNL^N) 

1  +  1x2. N£ .NX  +  ^)  GO  TO  302 
STAD.GT.C)  GO  TO  37 
ESTIM 

38 
S  E  S  T  1  M 
0. 

=  1^1  TO 

• 

0.   - 

S=  FMACH 
T=  0. 
=  0. 

r=  c.  <    ..      ■,  , 

(  IWPIT,600)  .  - 

( I  WRIT,  112)  V   ; 

x  +  1  .:.'.,.  Ki'  ^ 

1=  3,KX 

(IWRIT,61C)  A0(I),S0(I),S1(I),S2(I>,A1(I),A2(I),A3(I) 

( IwRIT, 600) 

( I  WRIT,  116) 

Y+2  .■  .  .r 

j=  3,MY  L  ■•'..:  .^.  - 

(IwRIT,6lC)  B0(  J)>B1(  J),B2(  J),B3(  J)  •  '  <''■■' 

5TA0.GT.0)  GO  TO  60  ■;•■:. 

ALF .bO.MHALF  )  MIT=  MIT*1.5 

.  i 

=LHALF*2**NHALF 

+  NY-KHALF 
3  -1 
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20 


Zl 


50 


51 


52 


I0« 

2  + 

KhALF 

11  ^ 

=  10 

+  1 

13  = 

2  + 

NX 

-KHALF 

12  = 

13  - 

-1 

INX=  13  -  10 
INY=  J3  -  Jl 
WklTE  (IWPITfoOO) 
WRITE  (IWRIT,12^) 
WKITE( lWRIT*b40) 
WRITEd  WRIT, 126) 
WkITE(IWPIT,6lO) 
WRITE( IWRIT>132 ) 
CALL  SFCONO(T) 
WRITEC IWRIT>66C) 
WRITE(IWRIT>128) 
NIT        =  r.  IT 
CALL  USTADl 
WRITE (IWRIT*650) 
IF(NIT.GE.MIT )  60 
IF (RG.GT.CaV)  GG 
CALL  SECONDd  ) 
WRITE  (IwRIT,660)  T 
IF(  NHALF.GE. MHALF ) 
NHALF=    NHALF  +  1 
M  I  T  =  M  I  T  /  ? 
NX=NX   +NX 
NY»    NY+   NY 
CALL  COORD {NL>N) 
IF(IXi+IX2.MF.NX+4) 
CALL  PEFIN 

Gu  ro  it 

USING  THE  STE  ADY 

WRITE  (iWRlT^bOO) 

WRITE  (IWRIT, 12^) 

WRITE  (  iWRIT^b'^O) 

WRITE  (lWRIT,12e) 

WRITE  (IWRIT^blO) 

CALL  SECaND(T) 

WRITE  (IWRIT,660) 

WRITE(IWRIT,129) 

MT=0 

ALT'U. 

PI         =  2./(l 

IF(  PlOl.  ECO.  )  GG 

Pl=  PIO 

IF  (  NHALF. EG.  1) 

IF(  NHALF.EQ.2. 

IF(NHALF.EC.3) 

P2=  P20 


I N  X  ,  I N  Y 

FMACH, AL 
DT 


+  1 


NIT*RG,IG>JG,KS 

TO  21 
TU  20 


CJ  TO  22 


GO  TC  302 


MODE  TO  GENERATE  THE  INITIAL  DATA 


NX,  NY 


FMACH, AL 


+  (2 

TO  51 


/PIO   -1.  )*.5**NHALF ) 


Pi=  PIO  +  PlOl 
)  Pl=  PIO  +  P102 
Pl=  PIO  +  P103 


P3 
Jl  = 
J3  = 
10  = 
13  = 


=  1 


3 

NY  +  2 
3 
NX  +  1 
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53 


5^ 


b5 


57 


58 


11  = 

12  = 
STEA 
NIT 
CALL 
WKIT 
IF  (N 
IF  (« 
CALL 
'wkIT 
IF( 
CALL 
CALL 
WRIT 
Wr^IT 
WRIT 
WRIT 
OU  5 
wRIT 
WRIT 
CALL 
WRIT 
CALL 
IROU 
CALL 
IPLG 
IF(  I 
GO  T 
NHAL 
i^iX=N 
NY  = 
CALL 
IF(  I 
CALL 
MIT  = 
GG  T 
STEA 
MT  = 
ALT  = 
ALTT 
FhAC 
FMAC 
CETA 
CETA 
CALL 
M1T  = 
WRIT 
WRIT 
WRIT 
KHAL 
J3=" 
J2  = 
10  = 


J3-1 
10  + 
13-1 

[JY  I 

STE 
E(IW 
IT. 6 
R.GT 

SEC 
E  (  I 
NHAL 

SVE 

FOR 
E  (  I 
E  (  I 
L  (  I  W 
t  (  I 
6  1  = 
E  (I 
E  (  I 

CPL 
E  (I 

SCH 
TE  = 

PSU 
T=  0 
ST  AD 
0  58 
F» 

X   + 
NY 

coa 

Xl+I 
SPE 

NIT/ 

G  1^ 

DY  I 

0 

0. 

=  0. 

hS  = 

HT  = 

T=  0 

TT  = 
EST 
MIT 

E  {  I  W 

E  (IW 

E(IW 

F=LH 

3  ♦ 

J3  - 

2  + 


TERATION  USING  THE  STEADY  MODE 

=  NIT   +1 
ADY 

RIT,()70)  NlT,RG*IG*JG*RP»IR*JK*TAIj*Pl*P2/P3*NS 
E.MIT)  GG  TO  5  5 
.CGV)  GO  TO  b^ 
aND(T) 
'.v  t^  1  T  ,  6  6  0  )  T 
F.LT.MHALF)  GO  TG  57 
LG 
CE 

WR IT>  60C ) 
W  R  I  T  ,  1  8  2  ) 

RIT>610)  AL*FMACri,CETA,CL,CD>CM 
WRIT» 18^  ) 
IXl, 1X2 

WRIT>610)  XP(I),YP(I),SV(I)»SM(I),CP(I) 
WPIT^  600) 
OT 

WRIT, 60C ) 
ART 
1 
PE 

I.GT.O)  GO  TLi  303 

NHALF  +  i 
NX 

+   NY 
RD(NL ,N) 

X2.Nt.NX  +  '^)  GG  TG  302 
FIN 
2 

TERATION  USING  UNSTEADY  MODE 


FMACH 
0. 

• 

0. 

IM 

0/2 

RIT,&00) 

RIT,13^) 

RIT,12^) 

ALFi'2**  NHALF 

NY-KHALF 

1 

KHALF 


'  1  0  <-  M  ' 
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59 


60 


22 


23 


35 


II  = 
13=  2 
12=  I 
INX" 
INY  = 
WRITE 

;rite 

WRITE 

WKITE 

CALL 

kvRITE 

WRITE 

NIT 

CALL 

WRITE 

IF(NI 

IF(  N 

IF  (KR 

CALL 

WRITE 

LNSTE 

N  I  T  =  0 

INITI 

NITS  = 

JSTEP 

CALL 

CALL 

WRITE 

WRITE 

WRITE 

WRITE 

DO  23 

V.PITE 

W  fi  I  T  E 

CALL 

WRITE 

CALL 

IROUT 

CALL 

IPLOT 

IF  (IS 

IROUT 

CALL 

IPJUT 

CALL 

CALL 

WRITE 

IF  (IS 

WRITE 

ISTEP 

WRITE 

MSTEP 

IJKLM 

K  F  a  P  C 


10  +1 

+  NX  -KhALF 
3  -1 
13  -  10 
J3  -  Jl 

(iWRlT^b'tO)  INX>INY 
(IWRIT*126) 
(IWRIT*610)  FMACHfAL 
(IWRIT,132)  DT 
SECONO(T  ) 
(IWRIT»6bO)  T 

( I  WRIT, 128) 

=  NIT   +1 
USTADI 

(IwRIT,650)  NIT,RG, IG,  JG,NS 
T.LT.3)  GO  TO  59 
IT.GE.MIT)  GO  TO  bO 
.GT.COV)  GO  TO  59 
SECCND(T  ) 

(IWRIT>56G)  T 
ADY  CALCULATION TIME  r^ARChlNG 


AL 

1 

=  0 

VELD 

FORCE 

(  IWR 

(  IWR 

(  IWRI 

(IWP 

I=IX 

(IWR 

(  IWR 

CPLOT 

(  IWR 

CHART 

E=  i 

PSUPE 

=  0 

YM.GT 

E=  IR 

LORD 

E=  IR 

SONIC 

SECGN 

(IWR 

TADI  . 

(6,60 

=  PER 

(IWRI 

=  1ST 

N=  MI 

=  UK 


DATA 


IT, 600) 

IT,  182) 

T,olO)  AL,FMACH,CETA,CL,CD,CM 

IT, 18^) 

1,  1X2 

IT,61G  ) 

IT, 600) 

IT, 600) 


XP( I ) ,YP(I),SV(I  ), SM(  I),CP (  I) 


.0.AND.ALS.EQ.0..AND.AMPLA.EC,0..AND.AMPLC.EO.0.) 
OUTE  +  1 

aUTE+  i 

D(T  ) 

IT, 660)  T 

GT.O)  GO  TO  303 

0) 

IGD*2.*PI/ (FREQR*DT) 

T,136)  ISTEP 

CP/  IPSURE 

N0( 800,  ISTEP) 

LMN/  IPSURE 


GDT035 


151 


2^ 
Zb 

Zb 
3^ 
33 
27 


28 


29 


36 

30 
31 

32 

39 

^1 

301 


JCHECK= 
IF (MOD  (M 
KFaRC=  K 
GU  TG  Z<* 
r!STEP=  I 
LSTEP:^  I 
KSTEP=  1 
KK3TEP= 
CALL  UST 
IF(PERID 
IF(MOD{K 
&D  TL  3  0 
IF(KSrEP 
GG  TC  S't 
IF (  IkOUT 
CALL  ROU 
IkaUTE= 
CALL  PLO 
IPLUT=  - 
Ir:GUTE  = 
CALL  VEL 
JSTEP=  J 
CALL  FOK 
WRITE  (I 
h'klTE  (  IW 
WRITE  (I 
WRI  IE  (Iw 
wkITt  (  I 
Ou  29  1= 
WRITE  (I 
CALL  PSU 
IPLUT»  0 
IF(  ISYM. 
IRQUTE= 
CALL  LOR 
IRGUTE= 
CALL  SON 
GO  ro  32 
IF{ MGD( K 
GO  TO  32 
CALL  VEL 
JSTEP=JS 
CALL  FOR 
KSTEP=Ki 
IFd^DDiK 
IF(KSTEP 
GO  TC  Zb 
IF(  ICHEC 
GO  TO  39 
CALL  TRA 
Ji=  3 
KHALF*  3 
J3=»  3  + 
10=  2  + 


lSTfcP*.2^ 

STEP,KFuPC) .EQ.O)  GO  TO  25 

FbkC-l 

STEP/ (KFORC*IPSURE  ) 
STEP    +    !) 

• 

FLOAT (  ISTEP) *.b 

ADi 

D.GE.2.)     GO    TG    3  3 

STEP^MSTEP)  .EQ.G)     GO    T[]    27 

.LT.KKSIEP)     GO    TO    30 

E.LT'.25)     GO     TG    2u 

TE (faLaUTPUT>2LLP ) 

0 

T(0.f0.^9  99) 

1 

IROUTE  +1 
u 

STEP  +  1 
CE 

WR IT, 600) 

RIT,13fc)  GTIM^KSTEP 
W  R  I  T ,  1 8  2  ) 

RIT^blO)  AL,FMACH,CETA,CL,CD, Cf 
wPIT,ie^) 

IXI,  1X2 
WRIT, 610)  XP (  I),YP( I ),SV( I),SM(  I),CP (  I) 

RE  .   :  -     ' 

GT  .0. AND.ALS.tQ.O. .AND.AMPLA.EC.O. . AND.AMPLC.EQ.O. )  GDT036 
IRJUTE  +  1 

D 

IRGUTE+  1  ■•-'•  ..:■.:,,.,■,■  ^ 

STEP,NSTEP).EO.O)  GO  TO  31  ''  '  '       '  '    '  "  ''  '   "   '  ■■''' 

0  .'•  ■' 

TtP  +  1 

CE 

TEP  +1  \'-  •    ^  '    ■■:     ■   ■'-■■'■ 

STEP, JCHECK)  .Eu.O)  GO  TO  ^1 
.GE.LSTEP)  GC  TO  301  '  •  ' 


K.GT.O)  CALL  CHEKPTX(VAR  ) 

CE 

*  2**MHALF 

NY-KHALF 

KHALF 


1  >''. 
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303 
302 
112 


13 
CA 
CA 
ST 

FU 


1 
2 
3 

116  FO 

1 

2 
12^  FD 
126  FQ 

128  FQ 
1 

129  FG 
1 

2 
3 
132  FG 


13A 
136 
138 
182 


FQ 
FQ 
FG 
FG 


=  2  + 
LL  GR 
LL  PL 

OP 
RMATC 


RMAT( 


RMAT( 
RMAT( 
RMAT( 
5H 
RMAT( 


RMAT( 
RMAT  ( 
RHAT( 
RMAT( 
RMATC 


NX    -KHALF 
ID 
GT {Q.f 0.>999  ) 


-^IHOMAPPED 
15hO 
15H 
1!?H 

IbHOY  STRE 
15H0 
15H 

15H0  HOii  D 
13 HO  MA 
lOHOITERAT 
>  lOhS 
lOHOITERAT 


COORDINATES  AND  X 


,lbH 

>  l3H 

) 
FACTORS/ 

,15h 

) 

IVISIGNS»15H 
C H  NO    ,lbU 


X 
YPP 

A3 
TCH 

Y 
B3 


STRETCH 
Y 
Al 


FACTORS/ 
/15H 


15H 


VER 
ANG 


Bl      >15H 

DIVISIONS ) 
OF  ATTACK) 


YP 
A2 


B2 


CORRECTION  »5H 


) 


16^  FGKMAK 

1 

2 

3 
500  FG 
510 
530 


600 
610 
620 
630 
6^0 
650 
660 
670 


FQ 
FO 
FG 
FG 
FG 
FG 
FO 
FO 
FO 
Fu 
EN 


RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
RMAT( 
PMAT( 
D 


lOH  CI 

lOH  SG 

1H0»*T 

1HG»*G 

1H0>*G 

IH  >*T 

15H0  A 

15h 

36hOCG 

26H  AN 

13H 

15H 

IX) 

bFlO.7 

20A^) 

IHl) 

8F13.A 

8E15.5 

1HC»2C 

18,711 

HO,  bi 

15H0Ca 

I10»Ei 


IGN,15H 

GNIC  PTS 

IGN,15H     CCRkECTICN 

15H      RESIDUAL 
ATN,1CH  REL  FCT  1,10H 
PTS  ) 


I  >5H 


RCUL 

NIC 

IME  STEP  =   ♦,F15.10) 

NSTE 

NSTE 

IME  = 

NG  G 


,5H 
,5H 
REL 


I 

I 

FCT 


,5H 
,5H 
2>10H 


J  , 
J  , 
REL  FCT  3, 


QRDI 
0    PR 


MA 


ADY  ITEkATIGN*) 
ADY  STEPS  =  *,  5X,I10) 
*,5X,F10.5,5X,*STEP  =  *,5X,I10,/W 
F  ATTACK, 15h   FLIGHT  SPEED  ,15H 

CL      >15h         CL      ,15H 
NATES  OF  INTERPOLATED  AIRFOIL* 
ESSURE  DISTPIBUTIGN/ 

X       >15H         Y       ,15H 
CH  NO    ,15ti  CP       ) 


FLIGHT  ANGLE  * 
CM       ) 


V/VO 


) 


) 
) 

A4) 
5) 

5.5, 
MPLT 
5.5, 


215,110) 

ING  TIME,F10. 3, lOH    SECONDS) 

215, E15. 5,2I5,4F10.5,I10) 


SUBROUTINE  COGRD(L,N) 

SETS  UP  MODIFIED  PARABJLIC  COORDINATE  SYSTEM 

COMMON/ A/  GM(132,36),G{132,36),GN(132,36),S0(132),S1(132),S2(132) 

1  ,A0(132),A1(132)»A2(132),A3(132),B0(36),B1(36),B2(36) 

2  ,B3(36),NX,NY,IX1, 1X2, KS YM, F MACH, AL PHA, C A, S A, F MACH2 

3  , AL,UTIM,CB,SB,NS,RG, IG,JG 

COMMON/C/  XP ( 260 ), YP ( 2b0 ), Dl ( 260 ),D2( 260), 03(260) 
CGMMGN/D/  SLGPT, TRAIL, SCAL 
COMMQN/F/  XR,YR,KS,XS{500),YS(500) 


153 


1?. 


23 


Z^ 


22 


CUf^M 
Cu.'iM 
CGflM 
PI 

DX  =  4 
CY 

OXY  = 
DYY  = 
DXX  = 
KX  = 
MX  = 
KY 
MY 

S=  C 
1=  C 
XT£  = 
SCAL 
DJ  1 
XO 
YO 
R 

ANGL 
IF  ( 
IF  ( 
R 

ANGL 
XP(  I 
YP  (  I 
DO  2 
XX  = 
b 


CN/H/ 

UN/I/ 
LN/0/ 

.  /N'X 

D  X  *  0  Y 
DY*OY 
DX*DX 
NX+1 
NX  +  2 


UDRS 
CGRT 
1. 


OX»OY,DT>DXX»OYY,DTT>OXY>DXTfDYT»TSR 

X (260) f Y(260) 

COaRS^COCRT 

=  3.1<.15q265358979 

=  1  ./NY 


=  NY 

=  t;Y 


+  1 

+  2 


=  . 50001*XTE**2/ (X(N)   -XR ) 


2  1  =  1, N 


I.LT. 
I  .GT  . 


L  . 
L. 


S 
S 

s 

c 

AN 
AN 

s 


) 

) 

2  1  = 
(  1-2) 


3f 


IF 

Sa 

CX 

XC 

XI 

X2 

Xi 

GO  T 

IF  ( 

A 

XO 

XI 

a2 

IF  ( 

IF  ( 

AC(  I 

Al  (I 

A2(  I 

A3(I 

IXl 

UO  3 

YY« 

B 


R 

R 

KX 

X 

=  1 

(  ABS(XX  )  .G 

•=  S 

=  C 

=  X 

=  1 

=■  S 
=  ( 


CAL*(X ( I )   -XR) 

CAL«( Y( I )   -YR ) 

QRT(XC*XO   ♦YO^YO) 

MPLA(XO*YO) 

D.  ANGL.LT.  .  !j*PI  )  ANGL 

D. ANGL.GT.l. 5*PI )  ANGL 

QRT(R   +R) 

5*ANGL 

*CCS( ANGL) 

*S  IN (ANGL) 


ANGL 

=  ANGL 


+  PI 
-PI 


+  PI 
-PI 


-2. 

• 

T.XTE  ) 
IN( PI* 
DS(PI* 
X  +S* 
./(i. 
*PI*SX 

1.  +s 


n  2  4 

XX. LT 


.0, 


) 

=  1 
=  b 
=  ( 

XO.LT.XP( 
XO.LE.XP( 
)  =  X 
)  =  . 
)  =  X 
)  =. 5*X2/0 
=  I 
3,  MY 
*UY 
=  1 


.   -(  ( 

*XTE 

1.   +S 

2. ♦(XX 

1)  )  IX 

N  )  )  IX 

0 

5*Xi/D; 

1*X1 

X 

XI    +1 


GO 
XX/X 
XX/X 
XTfc* 

+  s* 

*X1/ 
)*X1 

1. 
XX 

+  (XX 

)  *A* 

-d 

1  = 

2  = 


TO  23 

TE) 

TE  ) 

SX/ (PI*( 1. 

(1.   +CX)  ) 

XTE 


+  S  )) 


-B*XTE)/(2.   -XTE))**2  '   " 

-B*XTE )/(A*  (1.   +S)  )       ■  '  '    ^ 
A/ (2.   -A) 
♦xrF)*(^.   -A)/(A*(2.   -A)*(2.   -XTE)**2) 

I 


2  J  = 
(  J-3) 


-YY*YY 
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32 


^2 


52 


62 


72 


Yl 
*  (4 


GT. 
GT. 


F  ( 
L  ( 


Yl 

BO(J  ) 

Bl(  J) 

B2( J ) 'Yl* 

B3(J  )  =  -YY 

ANG 

ANGi 

IF  (ANGI. 

ANG2 

IF  (ANG2. 

ANGI 

ANG2 

Tl 

T2 

CALL  5PLI 

CALL  INTP 

XI 

S0{2)=  0. 

S0{f1X)=  0 

M 

A 

C 

DO 

XX 

XO 

YO 

R 

ANGL 

IF  ( ANGL. LT 

R 

ANGL 

SO(I 

M 


8*6/ ( (2. 
T*YY/B 
. 5+Yl/DY 


-B)*T) 


.-B)/ {B*(2.-B)*UY ) 

ATAN( SLOPT ) 

CMPLA(  (X(l)   -Xk)*  (  Yd  )   -YR)  ) 

PI)  ANGI  =  ANGI   -Pi   -PI 

CMPLA(  (X(N)   -XR)>  (Y  (N)   -YR)  ) 

PI)  ANG2  =  ANG2   -PI   -PI 

=  ANG   -.5*ANG1   +.5*TRAIL 

=  ANG   -.5*ANG2   -.5fTRAlL 

«  TAN(ANGl) 

=  TAN{ANG2) 


1,N»XP,YP>D1»D2»D3>1>T1,1,T2»0,0.,IND) 
IXl,IX2,A0,S0»i»N>XP*YP»Dl»u2#D3»0) 
X( 1)   -.75* (X(l)   -X(L)  ) 


XI   -1 

LL1PT*(X{1)   -XI) 
./  (X(l)   -XI) 


'^'d    1=  3»M 


+  XR 


5*A0( I )**2/SCAL 

CAL*(XX   -XR) 

CAL*{Y{i)   +A*ALOG(C*(XX   -XI)) 


-YR) 


) 


A 
C 

DG 

XX 

XG 

YO 

R 

ANGL 

IF  ( 

R 

ANGL 

S0{  I 

SCAL 

DG  6 

DSI 

OS  I  I 

SKI 

S2  (  I 

DC  7 

XP(I 

YP  (I 


52  I«  M. 


=  S 

~    • 

=  R 
=  I 
=  S 
=  1 
KX 

• 

S 

s 
s 
c 


URKXO^XO 
MPLA(XO,YO) 
5*PI)  ANGL 
QPT ( R   +R) 
5*ANGL 
*SIN( ANGL ) 
X2   +1 
LDPT*(X(N) 
./  (X (N)   -Xi  ) 


+YO*YO) 


=     ANGL       +PI        -t-PI 


■XI) 


+  XR 


5*A0(I )**2/SCAL 
CAL*(XX       -XR) 
CAL*(Y(N)        +A*ALGG(C*<XX 
QkT(XO*XU       +YO*YO) 
MPLA(XO»YO) 


-XI))        -YR) 


ANGL .GT.l .y*PI )     ANGL 
=     SQRT(R       +R) 


=    ANGL       -PI       -PI 


) 

2    1  = 

=  (S0( 

) 

)  =  A2  ( 

2    1=1 

) 

) 


3* 

1  + 

3 

I  ) 

XI 


K 

1 
KX 

S 
1) 

A 
*D 


5*ANGL 

*SIN{ANGL) 

./SCAL 

0(1+1)       -SO(I-l) 

-2.*S0(I)+30(I-1))/DXX    +A3(I)*DSI 

1(I)*0SI 

SII 

X2 

5*SCAL*(A0( I)**2       -S0(I)**2) 

CAL*AO( I)*SO( I )       +YR 


+  XR 


155 


KETUPN 
END 


SUBROUTINE   bPLIF(M,N,i,F>FP,FPP>FPPPfKM»VM»KN*VN*MODEfFQM*IND) 
C      SPLINE  FIT  -  JAMESON 

C       INTEGRAL  PLACED  IN  FPPP  IF  MODE  GkEATEF  THAN  0 
C      IND  SFT  TO  ZEPD  IF  DATA  ILLEGAL 

DIME^SIDN    S(1)^F(1)>FP(1),FPP(1),FPPP(1) 

IND         =0 

K  =  IABS(N  -A) 

IF  (K  -1)  ei^ei^i 

IK  =  (N   -M)/K 

I  =  M 

J  =  M   +K 

DS  =  S ( J)   -S  (  I) 

D  =  DS 

IF  (DS)  11»81>11 

il  OF  =  ( F( J  )   -F ( I ) )/DS 

IF  (KM   -2)  12,13,  lA 


12  U 
V 

Gu  TQ  25 

13  U 
V 

GC  TO  25 
i^  LI 

V 

GO  rC  25 
21  I 

J 


=  3.*(DF   -Vri)/DS 

=  0. 

=  VM 

=  -1. 

=  -DS*VM 

=  J 

=  J   +K 

=  S  (  J  )   -S  (  I) 


IF  (  C*DS  )  81,  81,23 

23  DF  =(F{J)   -F(I))/US 

6  =  l./(OS   +DS   ♦U) 

U  =  P*OS 

V  =  e*(b.*DF   -V) 
25  FP (1  )  =  U 

FPP( I )  =  V 

U  =  (2.   -U)*DS 

V  =  e.*DF   +DS*V 
IF  (J  -N)  21,31,21 

31  IF  (KN   -2)  32,33,3^ 

32  V  =  (6.*VN   -V)/U 
GO  TO  35 

33  V  =  VN 
GO  TQ  35 

3^  V  =  (  DS*VN   +FPP(I  )  )/(l. 

35  B  =  V 

0  =  OS 

^1  DS  =  S (J)   -S(  I) 

U  =  FPP(  I  )   -FPd  )*V 


+  FP(  I) ) 


156 


31 


61 


71 


81 


(J       -M  ) 


FPPP  (  I  ) 

FPP(  I  ) 

FP(  I) 

V 

J 

1 

IF 

I 

FPPP(N) 

FPP(N) 

FP(N  ) 

IND 

IF    (MODE) 

FPPP (J  ) 

V 

I 

J 

DS 

U 

FPPP( J) 

V 

IF     (J 

RfcTUt'N 

END 


-F{I))/DS       -DS*(V       +U      +U)/b. 


-N) 


=  (V       -U)/DS 

=  U 

=  (F(  J) 

=  U 

=  I 

=  I       -K 

-tl*  bl*  -tl 

=  N        -K 

=  FPPPd  ) 

=  B 

=  DF       +D*(FPP( I ) 

=  1 

81,81>61 

=  FQM 

=  F?P( J ) 

=  J 

=  J        +K 

=  S  (  J)       -S (I) 

=  FPP(J  ) 

=  FPPP( I )       +.5*DS*( F( I) 

=  U 

7i»eiWl 


+B       +B)/6. 


+F(J)       -OS*DS*(U       +V)/12.) 


3 
ti 

11 

12 
13 

21 
31 


FUNCTION  CMPLA(X^Y) 

ANGLE    OF    CQf'.PLEX    NUMBER     X     +I*Y 
PI  =    3.1413y26i»3f)8979 

IF     ( ABS(Y  )       -ARS(X)  )     1,1*11 


IN  RANGE  a,  TU  2,*PI 


SHIFT 
IF  (X) 
SHIFT 
IF  (Y) 
SHIFT 
C  M  P  L  A 
GC  TO 
SHIFT 
IF  (  Y) 
SHIFT 
C  i^i  P  L  A 
GO  TO 
C^5PLA 
RETURN 
END 


=  PI 
^,21,2 

=  0. 
3>^,<. 

=  2.*PI 

=  SHIFT 
31 

=  .5«PI 
12,12*  13 

=  1  .i>*PI 

=  SHIFT 


+ATAU ( Y/X ) 


-ATAN (X/Y ) 


31 


=  0. 


C 

c 

c 


SUBROUTINE   INTPL(M,NI,SI,FI,M,N,S,F,FP>FPP,FPPP,MODE) 

INTERPOLATION  USING  TAYLGk  SIRIES  -  JAMESON 

ADOS  CORRECTION  FUR  PIECEWISE  CONSTANT  FOURTH  DERIVIATIVE 

IF  MODE  GREATER  THAN  0 

DIMENSION    SI(1),FI(1),S(1),F(1),FP(1),FPP(1),FPPP{1) 
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11 

13 

15 
21 


23 

31 

33 

35 
37 


K  =  IABS(N   -M) 

K  =  (N   -M)/K 

I  =  ,^ 
MIN         =  MI 
NIN         =  NI 

D  =  S (N)   -S ( M) 

IF  (D*(Sl(iMI)   -SI(MI)))  11*13^13 

MIN         =  NI 

iM  I  N         =  i1 1 

KI  =  IABS(NIN   -MIN) 

IF  (KI)  21.21»15 

KI  =  (NIN   -MIN)/KI 

II  =•  .1IN   -KI 
C           =0. 

IF  (MUDE)  il>il^23 
C  =  1  . 

=  I  I   +KI 
=  SKID 
=  I   +K 
(1   -N)  3b>37,35 
(D*( S( I)   -SS) )  33>33,37 
=  I 

=  I    -K 
=  SS   -S( I ) 

=  C*(FPPP(J)   -FPPP( I)  )/ (S  (J)   -S(I)) 
=  QUARP(SS*F(I)>FP(I)>FPP(I)»FPPP(I),FPPPP) 
-NIN)  31»A1>31 


'il 


II 

SS 

I 

IF 

IF 

J 

I 

SS 

FPPPP 

FI  (II) 

IF  (II 

kETURN 

ENO 


FUNCTION 

EVALUATES 

QUARP 

UUARP 

UUARP 

QUARP 

RETUPN 

END 


QUAPP(OS^F*FP,FPP,FPPP, FPPPP) 
FIRST  FOUR  TERMS  QF  TAYLOR  SERIES 
=  FPPP   +,25*DS*FPPPP 
=  FPP   +DS*UUARP/3. 
=  FP   +. 5*DS*UUARP 
=  F   +0S* QUARP  .  ■  '  , 


-  JAMESON 


-  .'•'  ,1 


SUBROU 

PARAQQ 

SO 

FPO 

FPP 

FPP 

FP 

F 

RETURN 

END 


TINE   PARAF(S1,52,S3,F1,F2^F3,F^FP»FPP) 
Lie  FIT  -  JAMESON 

=  .5*(S3   +S1) 

=  (F3   -F1)/(S3   -SI) 

=  (F3   -F2)/(S3   -S2) 

=  2.*FPP/ ( S3   -SI  ) 

=  FPO   -FPP*SO 

=  F2   -S2*(FP0   +FPP*(.5*S2   -SO) 


-(F2   -F1)/(S2   -SI) 
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SUBRDUTINE  SESTIM 
C       STEADY  ROUTINE 
C      INITIAL  ESTIMATE  OF  REDUCED  POTENTIAL 

COMMON/ A/    GM(132^3o)>Gll32*36),GN(132,36)*S0(132)>Sl(132),S2(132) 

1  ,A0(132)*AI(132)>A2(132)»A3(132)>B0(36)/Bi(36)*B2(36) 

2  ,B3(36)»NX,NY,IX1,IX2>KSYM,FMACH,ALPHA,CA,SA,PMACH2 

3  ,ALf UTIM,CB>SB>NS*RG>I&^JG 
CUMMON/M/     P1>P2>P3>TAU 
COMMON/WAKE/     NIT,WG(132) 

IX=  NX+^ 
IY=  NY+^ 
KX=     NX+1 

MY  =  NY   +2 

Cb=  CaS(ALPHA) 
S6=  SIN(ALPHA) 
CA=  FMACH*CB 
SA=  FMACH*S8 
TAU=  0. 
DO  12  I»  1,   IX 
DO  12  J=  l^IY 
12  G(  I»  J)     =0. 
DO  Z2    I=IX1»IX? 

uo  = 

8IS=    Bl(3)*(l. 

22  G{ I»2)=   G( I, ^)-(CA*SO(  I  ) 
DO  23  1=  IX2»KX 
M=  NX+^  -I 

23  WG(I )=  G (I*3)-G(M,3) 
RkTUKN 
END 


CA*AG(I)   +SA*SO(I) 
+  S1(  I) **2) 

-SA*AO( I )   +U0*S1(I) ) /BIS 


SUBROUTINE  SREFIN 
C      STEADY  ROUTINE 
C      HAL\/bS  MESH  SIZE 

COMMON /A/  GM(  132^30  )>G(  132,  36  ),GN(  132*36), SO  (132), SKI  32  ),S2(132) 

1  ,A0(132),A1(132),A2(132),A3(132),B0(36),B1(36),B2(36) 

2  ,B3(36),NX>NY,IX1,IX2,KSYM,FMACH,ALPHA,CA,SA*FMACH2 

3  , AL,UTIM,CB,S6,NS,RG,1G,  JG 
COMMON/WAKE/  NIT,WG(132) 

KX=  NX+1 


MX  = 

KY 

MY  = 

IY» 

LX« 

LY  = 

DO 


NX +  2 


NY  +  3 

NX/ 2 

NY/2 

22  K  = 


=  NY   +1 
NY  +  2 

+  2 
+  3 
2,LX 


1=  LX+2-K 
11==  (I-2)*2  +2 
DO  22  KK=  3,LY 
J=  LY+3-KK 
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+G( I,J-1) ) 
+  G(I-1,J  )  ) 


JJ*  ( J-3  )*Z    +3 
22  G(II/JJ)=       G(I/J) 

DO  42  1=  2>MX,2 

DG  42  J=  4>MYi2 
42  G( I>J )     =  .5*(&( I, J+1) 

DD  3  2  J=  3»  lY 

u'J     3^  1=  3/KX»2 

32  G(  I»J  )     ■=  .  5*(G(I+1*  J  ) 
DO  33  1=  IX2^KX 

ri=  NX +  4- I 

33  wG(I)=  G(I»3)-o(M>3) 
DU  b2  1=  IX1>  1X2 

GI=  G( I+l>3  )-G(  I-l»  3) 

U0=        Al(i;*GI  +CA*AO(I)   +SA*SO(I) 

e.IS=    3i(3)*(l.    +S1(I)**2) 
52  G(I»2)=  G(I>4)     -(CA»SO(I)  -SA*AO(I)  +U0* S 1 (  I  )  )  /  B  IS 

N  =  1X1   -1 

DO  62  1=  3fN 

iM  =  N  X  +  4  -  I 
fc2  G(M*2)=  6(1*4)  +  WG(K) 

N  =  1X2   +1 

DU  6^  1=  N*KX 

M=  NX+4  -I 
64  G( M,2)=  G(  I>4  )   -WG(  1  ) 

RETURN 

END 


SUB 
STE 
CAL 
CUM 

1 

2 

3 
COM 
COM 

1  *A 

2  ,F 
LQM 
AAO 
DG 
Y=B 
H  =  S 
GI 

G  J  = 
U 
V 
00  = 

u  =  s 

IF 
SV( 


ROUTINE  SVfcLU 
ADY  ROUTINE 


1  .  +  -_ 
12  1=1X1,1X2 
0{3)+S0(  I) 
QPT(AO (  I)**2  +SG(  I  )**2) 

G(I+1,3)-G( 1-1,3) 

G( I,4)-G(  1,2) 


U*U+V*V 
QPT(Q3) 
(U.LT .0. ) 
I) 


-G(  1,2) 

=(A1(I)*  GI   -S1{I)*B1(3)*  GJ   +C 

=  (31(3)*  GJ   +SA*AO(I)   -CA*Y)/H 


+CA*AO(I)   +SA*Y)/H 


=  -& 
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12 


AA«  AAO  -.2*QQ 

AA=  A6S(AA) 

A=  SORT(AA) 

SM(i  )«  Q/A 

CP{ I)=    ( AA**3.5 

RETURN 

ENO 


-1. ) /( .7*FMACH2) 


C 
C 


2 
11 


14 


12 


610 


SUBRUUTI 
STEADY  R 
GENERATE 
CGMMQN/A 

1 

2 

3 
COMMDN/J 

1  ,AMPLM> 

2  ffREQR, 
CGMMCN/K 
DIMENSIO 
AAO  =  1. 
IwRIT 

K 

IF  (NY.G 

WRITE  (I 

FCjRMATC  i 

DC  12  1= 

J=J1 

N=  0 

N=  N+1 

Y 

HH=AO(I  ) 

H=SQRT(H 

GI=G{  I  +  l 

GJ=G(I»J 

U 

V 

aQ=u«u+v 

AA=  AAO 
AA=  A3S{ 
QA=   QQ/ 
IND{N)  = 
IF  (U.LT 
J=  J+K 
IF(J.LE. 
^.RITE  (I 
RETURN 
Fl1RMAT(1 
END 


NE  SCH 

OUTINE 

S  MACH 

/  GM(  I 

>AO(i 

f  B3(  3 

>AL^U 

/  RAD» 

FkEQPM 

IPSURE 

/  10,1 

N  I N  D  ( 

+  .2 

=  6 

=  NY 

T.32+K 

WRIT, 2 

^HOf^AC 

10,13 


ART 

NU  CHART 
32,3b),o(132,36),GN(132,36),S0(132),Sl(132)>S2(132) 
32),Al(lj2),A2(132),A3(132),80(36)#Bl(3fa),B2(36) 
b),NX,NY,IXl,IX2,KSYM,FMACH,  ALPHA*CA>SA*FMACH2 
TIM,CB,Sf^,NS*RG,IG,JG 

PI*ALj>ALT>ALTT,AMPLA*FREQRA»FASAGA,FMACHS,FMACHT 
,FASAGK,CETAS,CETAT,CETATT,AMPLC,FREQRC,FASAGC*CETA 

1,12,  13, Jl, J2,J3 

150) 

♦  FMACH2 

/32 

)   K  =   K    +1 

) 

H  r:0  CHART) 


=  S0( I )   +B0( J ) 

*A0(I  )+Y*Y 

H) 

,J)-G(I-1,J) 

+1 )-G(I,J-l) 

=(A1(I)*  GI   -SI ( I)*B1( J)*  GJ   +CA*A0(I) 
=  (BKJ)*  GJ   +SA*A0(I)   -CA*Y)/H 

*V 

-.2*QQ 

AA  ) 

AA 

iOO.*SGPT ( QA) 

.0. )  IND(N)  =  -IND(N) 

J2)  GO  TO  1-^ 

WRIT, 610)  (IND( J), J=1,N) 

X,32I4) 


+SA*Y)/H 
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SUB 
STE 
STc 

PAR 
G  I 
AND 
CUM 


1 
2 
3 


21 


36 


COM 
CUM 
CCM 
COM 
DIM 
J^  = 
I^  = 
11  = 

no 

113 

JU  = 

EE  = 

NS 

AAO 

kK  = 

IP 

OR 

i<G  = 

IG 

JG 

Rt=» 

IE  = 

JE  = 

Jl 

U2 

Cd 

D(  I 

J  = 

DO 

FX  = 

Y=S 

Hh  = 

M=S 

DH  = 

GI  = 

GJ  » 

U 

V 

AV 

AU  = 

S 

IE 

T 

IF( 


RGUT 
ADY 
ADY 
A3  0L 
S  TH 

I  5 
MON/ 


MGN/ 
MGN/ 
MCN/ 
MdiM/ 
ENSI 

J3 

13 

10- 
=  10 
=  13 

1./ 

1./ 


INE  ST 

RCUTIN 
T><ANSG 
IC  COO 
E  VELQ 
THE  RE 
A/  GM( 
^  A0( 
,33( 
*  AL» 
H/  UX» 
K/  I0> 
M/  Pl> 
STADI/ 
GN  C(l 
+  1 
+  1 
1 
+  2 
-2 
DXX 
DYY 

=  0 


EADY 

E 

NIC  PLTENTIAL  FLOW  EQUATIGiN  IN  SHEARED 

ROINATES  SYSTEM  SGLVED  BY  ROW  RELAXATION 

CITY  POTt.'sTUL  IN  THE  ABSOLUTE  FRAME 

DUCED  POTENTIAL  IN  THE  UNIFORM  MOVING  FRAME 

132>36)>o(i32,36)^GN{132»36),S0(132)>Sl(132)»S2{132) 

132),Al(132)>A2(13  2),A3(132),80(36)»Bl(36),B2(3o) 

3t)>NX,NY, IX1,IX2>KSYM>FMACH* ALPHA>CA*SA>FMACH2 

UTIMiCB/5DjNS>RG*IG»JG 

OY>0T*DXX,OYY>OTT,CXY,DXT>DYT>PDT 

I1^I2» 13, Jl, J2, J3 

P2>P3>TAU 

RK*  IR»  Jr(,  IRSTAD 
32)iD(132) 


=  1.  +  .2  *  FhACH2 


0 


0. 


0. 

0. 

0. 


I  )  = 
I  )  = 

J3 

32  I 
l.+S 
0(1) 
A0(I 
OPT  ( 
1./ 
G(I  + 
G(I» 


U  + 

(U.L 

AV. 


=  0 

=  c 

=  0 


=  2  .  /  P  1 
=  i  ./P2 


0. 
0. 

=  I 
1(1 
+  B0 
)»A 
HH) 
H 

i*  J 
J+1 


V* 

T.O 

1. 

LT. 


0,  13  ' 

)**2 

(J  ) 

0(  I  )  +  Y*Y 


)  -G(I-1>  J  ) 

)-G(If J-1) 

=(A1(I)*  GI   -SI ( I) *Pi( J)*  GJ   +CA«AO(I) 

=  (BKJ)*  GJ   +SA*AO(I)   -CA*Y)*DH 

=  V   -U*S1(I ) 

SKI)  .   •  .  .  ■ 

,  1.  ■     ->'*  >   '^■.^-  ■■■'  '     ''  -       '  ■  - 

. )  s  =  -1.  '  *'  '••  '■ 


+SA*Y)*DH 


0.  ) 


-1. 
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UU=U*U 

VV>=V*V 

UQ»UU+VV 

AA=  AAO  -.2*QQ 

A6  =  A1(I  )*B1( J  ) 

GII=(   G(  I  +  l» J)-G(I> J)-G( I> J) +G(  1-1, J ) )*Dt    +A3(I)*GI 

GIJ   =G(I+1*J+1)  -G(I-1*J+1)  -G(I+1*J-1)  +G(I-1»J-1) 

GJJ=    (G(i>J+l)  -G(IfJ)  -G(I*J)  +G(I,J-I))*EE   ♦  fl3(J)*GJ 

R=  -(AA   -UU) *S2( I)*bl( J)*  GJ 

1  +CA* ( VV-UU)-2.*UV*SA 

2  +  CC*(U*AO( I ) +V*Y )*DH 
IFCQQ.GE.AA)  GO  TO  33 

AXX         =  ( AA   -UU)*A2(I) 
AXY=-( AA*S1(I)+U*AV)*2.*A3 
AYY=(AA*FX-AV*AV)«B2(J) 
R=AXX*GII+AXY*GIJ+AYY*GJJ+R 


33 


Ai=    -2.*AXX*L'D 
131=     AXX*DD 
Cl=     AXX*DD 
YI  =  -R 
GO    TC    35 
NS 


66 
67 


6^ 
hi.' 


-Q1*AYY*EE 


NS 
S 
I 
IM 


+  1 


K 

IM  =     I        -K 

IMM  «     IM       -K 

L=  T 

JK=     J-L 

Ji^M=     JM-L 

AQ  =     AA/QQ 

BXX='VV*A2{  I) 

BXY=    -2.*AB*V*AU 

BYY=     AU*AU*B2(J) 

GNN=BXX*GIl+BXYtGIJ+BYY*GJJ 

IK     1MM.LT.2.0R. IMM.GT.  I^)     GO    TO    66 

GIIM=  (G(1>J)     -G(IN.iJ)     -G(IM^J)     +G  {  IMM>  J  )  )  *D0 

GU    TU    67 

GIIM=    Gil 

GIJM=  G(I,J)       -G(II'.,J)  -G(IfJM)  +G(IM,JM) 

iF(     JMM.LT.2.  Uk.  JM!1.GT.  J^)     GO    TO    bH 

GJJM  = 

GU    TO    6!? 

GJJM  = 

AXX 


+A3{ I)*GI 


(Gd^J)     -G(I^Jfl)    -G(I»JM)     +G(I,JliM)        )*td       +    B3(J)*GJ 


GJJ 

=     UU*A2(I) 
AXY*8.*S*T*U*AV*AB 
AYY  =  AV*AV*B2( J  ) 
GSS=AXX»GIIf1+AXY*GIJM+AYY*GJJM 


bt3  = 

Al  = 

Bl" 
Ci  = 
YI  =  -R 


*     (AC       -l.)»GSS       +AQ*GNN       ♦« 
(AQ-1, )*( AXX*DC     +     .5*AXY) 
AC*{-C2*BYY*EE       -2.*BXX*DD) 

+(A&-1.)*    2.*     (AXX*DD       +    AYY*EE       +AXY) 
AQ*3XX*0D-(1 .+S)*BH 
AQ*BXX*DJ       -(l.-Sl^BB 


37 


32 


^3 

4? 


32 
62 

32 
91 
92 
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lt-(AbS(R).LE.Rl^)  GG  10  37 

RR=ABS{ R  ) 

1  r'.  =  I 

JP=J 

A=  1. /(  AI-BI*C ( I-l  )  ) 

C ( I )=  CI*A 

D{  I)=  (Yl-3I=l--0(  I-i)  )*A 

CONTINUE 

1=  13 

CG  =  0  . 

Lj  '^2  y\=    I0fl3 

CG         =0(1)   -C( I) *CG 

IF  (AbS(CG).LE.A3  5(KG))  GG  TO  ^3 

kG= AcS ( CG) 

IG=I 

JG=J 

G(I>J)=    G(I^J)  +CG 

I  =1-1 

J  =  J-1 

1F(  j.GE.Ji)  GG  TD  21 

TTAU=  b( Ia2^3)-G( IX1>3) 

IF(  kSYiM.LE.J)  1AU=  TAU+  P  3*  (  TT  AU-T  AU  ) 

Dij  52  1=1X1,  1X2 

Gl=  G(I  +  i,3)-G(I-l>3) 

UU  =  AKD*  Gl   +CA*AO(I)   +SA*SO(I) 

BIS=    Dl(  3)*(1.  +S1(1)**?) 


-1 


+  TAU 
+  1 


G(  I>2)=  G(I»^  ) 

CCNTINUE 

N  =  IXl 

DG  6  2  1=  10,, M 

M=  NX+  ^  -I 

G ( M  ,  2  )  =  G (I  ,  ^  ) 

K  =     IXZ 

iJU  Iff  I  =  N,  13 

M=  iMX+  4  -I 

G{ M,2)=  &(  I>^ )-TAU 

IF(  FhACM.LT.l.  )  GO  TO  91 

DQ  dZ    J=  JlfJ^ 

3.*(G(IC>J)-6(I1,J)) 
3.*( G( 13, J )-G(  I2>  J)  ) 


(CA*50(I)-SA*AO(I)+UO*S1(I))/BIS 


G(  II,J  )  = 
G(  I^, J)  = 
i^E  TURN 
U  Lj  9  2  J  = 
G(II,J  )  = 
G(  K>  J)  = 
G(  II,J<t  ) 


G(  no,  J  ) 

G(  113, J) 


Jl^  J3 
-.b*TAL 
.  5*T  AU 
-.2b*IAU 


&(1^,J^)=  .25-*TAU 

KETURN 
END 


5U3RQUTINE   E  S  T  1 1"', 

INITIAL  ESTIMATE  GF  PUTENTIAL 
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Cj'IMCN/A/     GM(  13^,36),G(132,36),GN(132,36),S0(132),S1(132  ),S2(132) 

1  ,A0( 132 )>A1(132)>A2(132)»A3{ 132), B0(36), 61(36), BZ( 36) 

2  ,B3(36)»NX,NY,IXl,IX2,KSYr^,FMACH,ALPHA,CA,SA,FMACh2 

3  , AL, UTIM>CB,Sb,NS,kG, IG, JG 

COMMON /H/     DX,OY,DT,OXX,DYY,DTT,DXY,OXT,DYT,TSR 

CaMMUN/J/     kAD,PI,ALS,ALT>ALTT,AMKLA,FREOKA,FASAGA,Fr'.AChS,FMACHT 

1  ,AMPLM,FREQkM, KASAGM,CETAS,CETAT,CETATT,AMPLC,FREGPC,FASAGC,CETA 

2  ,FREQR, IPSUkE 
CQMMON/STADI/  RP , IR , J R , I RS T AD 
CQMMUN/WAKE/  NIT,WG(li2) 

KX=  NX+1 

MY  =  NY   +2 

CB=  COS(ALPHA) 

Se=  SIN(ALPHA) 

CA=  FMACH*CB 

SA=  FMACH*S3 

If  (IPSTAD.GT.O)  GO  TC  11 

DG  12  1=  3,KX 

DO  12  J=  3, MY 

&M(I,J)=  0. 

6N(I,J  )  =  0. 

G{ I, J)=  0. 
12  CLMTINUE 
11  JQ  22  1=  IXl,  1X2 

Y=SO(I )+80{3) 

X=  AU(I) 

Hh=  X*X  +  Y*Y 

DHH=  l./HH 

XT=  -.5*Y*( ALT+CETAT)  +  DHH*(CA*X  ♦  bA*Y) 

YT=  .b*X*  (ALT  +  CETAT  )  -UHH* ( C A ♦ Y- SA* X  ) 

VbN  =  hH* (XT*S1 (  I)  -YT) 

oI=  G( I+l,3)-G( 1-1,3) 

FX  =  1.  +  S1  (I  )**Z 

blS=FX*Bi(3) 

&XSXVB=  Al(  I)*GI*S1 ( I)  +  VBN 

G(I*2)  =  G(I,4)  -GXSXVB/BIS 

G(I,1)=  G(I,5)  -2.*GXSXVd/BIS 

22  CONTINUE 

DC  23  1=  1X2, KX 
K=  NX+^  -I 

23  WG(I  )=  G(I,3)-G(M,3) 
HKIN=  10. 

DO  13  1=  3,KX 
DO  13  J=  3, MY 
Y=SO(  I)  +  BO( J) 
HH=Y*Y+AO( I)*AO( I ) 
H=SUkT (HH) 
HX=  .5*H/A1(I) 
HY=  .i>*H/ei(J) 
riMl:M=  AMIN1(HMIN,HX,HY) 
13  CQNTINUE 

DT=  HMIN*TSP 

IUT=  2.*PI/  (FREQR*OT  )  <■  1 . 

IDT=  IDT/IPSUPt  +  1. 
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DT«  ^.*PI/ (IDT*I?SURE*FkEJR ) 

L)TT=DT*DT 

DXT=LX*DT 

OYT=DY*DT 

R  truRN 

END 


SUBROUTINE   REFIN 
HALVES  MESH  SIZE 
COMMON/A/  GM(132> 

1  , AG (132) 

2  ,33(3d), 

3  ^  A  L  #  U  T  I  h 
CUIMCN/H/  DX,DY,D 
COMMUN/J/  yAOjPIf 

1  >  AMPLM>FPEQRf1,F  A 

2  ,FREQRf IPSURE 
COMMCN/k^AKE/  NIT* 
DCT=  DT 

HMIN=  10. 

KX=  NX+1 

MX«  NX+2 

KY  =  NY   +1 

MY=  NY +2 

DO  13  1=  3»KX 

Ou  13  J=  3*.",Y 

Y=SO(  I)+EO(J ) 

HH  =  Y*Y  +  AO(  I  )*A0( I  ) 

H=SQkT(HH) 

HX=  .5*H/A1(I) 

hY=  .i)*H/bi(J) 

HMIN=  AMIM  (HMIN*HX,HY) 
13  CONTINUE 

0T=  hMIN*TSk 

IDT=  2.*PI/  (F 

IDT=  IDT/IPSU 

DT=  2.*PI/{I0 

DT=  AMIN1(DT> 
1^  RATIC=  DT/DDT 

UTT=DT*DT 

DXT=DX*DT 

DYT=DY*DT 

IY=  NY+3 

LX=  NX/2  +  2 

LY=  NY/2  +3 

DO  22  K=  2,LX 

1=  LX+2-K 

11=  ( I-2)*2  +2 

DO  22  KK=  3*LY 

J=  LY+3-KK 

JJ=  (J-3)*2  +3 


3b)/G(132»36)K,N(  132*  36),  50(132  )*S1(  132)^52(132) 

*A1(132)>A2(132)*A3(132)*B0(36),B1(36),B2(36) 

NX,NY,IX1,IX2*KSYM*FMACH*ALPHA,CA*SA,FMACH2 

,CB*5fc,N5>R&>IG*J6 

T*DXX,DYY*DTT>DXY,OXT,DYT,TSR 

ALS^ALT, ALTT*AMPLA,FREQRA,FASAGA>FMACHS,FMACHT 

SAGM,CETAS*CETAT,CETATT,AMPLC,FREQRC#FASAGC,CETA 

WG{132  ) 


RECR*DT)    +     1. 
RE    +     1. 

T*IPSURE*FREQR) 
DOT) 
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22 


&M(IWJJ)=  GM(I»J)*RATia 
G(II* JJ  )=       G(  I, J  ) 


42 


32 


33 


2*MX,2 
4*MY>  2 
.5*( GM( I>J+1)  + 

=  .5*(G(1> J  +  1  ) 
3*IY 
3>KX,2 
.5*(GM(I+1*J)  + 

=  .5*(G(I+1*J ) 
IX2>KX 


GM(I>J-1) ) 
+  G(  I* J-1)  ) 


GM( I-lf J) ) 
+  G(  I-lf  J  )  ) 


OD  42  I- 

DO  42  J= 

GM(I*J  )  = 

G( I> J) 

DO  32  J= 

DO  32  1= 

GM(  I> J)  = 

G(I»J  ) 

DO  33  1= 

M=  NX+4-I 

WG(I )»  G(I»3)-G(M>3) 

DO  52  I=IXlflx2 

Y=SO( I)+B0(3) 

X=  A0( I) 

HH=  X*X  +  Y*Y 

DHH=  l./HH 

XT=  -.5*Y* (ALT+CETAT)  +  OHH*(CA*X  +  SA*Y) 

YT=  .5*X*{ALT+CETAT)  -DHH* ( C A*Y-S A*X  ) 

VBN=HH*(XT*S1  (I  )  -YT) 

GI=  G( I  +  l>3)-G(  1-1^3) 

FX=1.+S1(I  )**2 

BIS=  FX*t31(3) 

GX5XV8=  Ai(  I) «GI*S1(  I)  +  VBN 

G(I*2)  =  G(I>4)  -GXSXVd/aiS 

G(I»1)=  G(I>5)  -2.*GXSxVB/bIS 

GM(1»2)=  GM(I»4) 
52  CONTINUE 

N  =  IXI 

OD  62  1=  3*iM 

M=  NX+4  -I 

GM(M,2)=  GM(1*4) 
62  G(M>2)=  G(I>4)  + 

N  =  1X2 

DO  64  1=  N*KX 

M=  NX+4  -I 

GM(M*2)=  GM(I»4) 
64  G(M*2  )=  G( I ^4  )   • 

RETUKN 

END 


-1 


W  G  (  f  ) 
+  1 


WG(I) 


SUBROUTINE  VELD 

CALCULATE^  SURFACE  VELOCITY  c,oo, 

COMMON/ A/    GM(132,36)>G(  132,36),GN(132,36),S0(i32),Sl(132),S2(132) 

1  ,A0(132)>A1(132)»A2(132)»A3(132)*B0(36)*B1(36)»B2(36) 

2  ,B3(  36)  ,NX,NY,IXl,IX2,KSYM,F.iACH,ALPHA,CA,SA,FMACH2 

3  , AL*UTiM,CB,S3>NS^RG> IGi JG 
COMMON /B/     SV(132)*SM(132),CP{132) 

C0.1M0N/H/     DX,DY>DT,DXX,DYY,DTT,DXY,OXTfaYT,TSR 

COMMON/ J/     RAD*PI>  ALS*ALT>ALTT,AMPLA,Ff<EQRA*FASAGA#FMACHS*FMACHT 
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12 


1  y 
I    * 

LU 
DU 
mA 

JG 
Y  = 
X  = 
HH 
DH 
XT 
YT 
H  = 
DH 
GI 
GJ 
GX 
GY 
U  = 

v  = 

QQ 
Uk 
VK 
GQ 
OK 
IF 
SV 
CH 
Fl 
Aa 
AA 
A  = 
SM 
CP 
RE 
FN 


Art  PL 
FKEU 
lii^DN 

r=  1 

0  = 

12 
80(3 

A0( 
=  X* 
H=  1 

=  .5 
SCR 
=  1. 
=  G( 
=  G( 


^,FP 
R>  IF 
/K/ 
./DT 
1. 

I  =  IX 
)+S0 
I) 

X  + 
./HH 
5*Y* 
*X*  ( 
T  (Hh 
/H 

I  +  l^ 
1,^) 


=  b 

GX 
GY 
=  U* 
=  X 
=  Y 
R  = 
=  S 
(Ufr 
(I  ) 
AIN 
T  = 
=  A 
=  A 
So 
(I  ) 
(I  ) 
TUR 
D 


1(3)* 
*DH 
*0H 
U+V*V 
T*H  + 
T*H  + 
U  R  *  U  R 
URT  (Q 
.LT.O 
=  GR 
=  XT* 
GM(  I, 
AO  -. 
BS  (AA 
RT  (  AA 
=  OR/ 
( 
N 


EQf'M,F  ASAGK,CEI  AS»  C  E  T  A  T  ,  C  E  T  ATT  »  AMP  LC»  F  REQRC  ,  F  AS  AGC  >  C  ETA 

SURE 

I0»  I1>I2*I3>J1*  J2>  J3 


1^  1X2 
(  I  ) 

Y*Y 

(ALT+CETAT)     +     DHH*(CA*X    +     SA*Y) 
ALT+CETAT)     -DHH* ( C A* Y- SA*X  ) 
) 

3  )-G(  \-\>  3) 

-G( 1*2) 

=    Aid)*    GI       -Sl(  I)*B1  (3)*    GJ 

GJ 


U 

V 

+     VP*VR 
QK) 
.)     GR=    -GR 

GX     +    YT*3Y 

3)     *DUT     +     CHAIN 

2*0G-. A*FIT 

) 

) 

A 

AA-!=*3.i       -1.  )/ (  .7*FhACH2) 


SUBROUTINE    FORCE  -'         ■      '  ■  '      '  ' 

CALCULATES    FORCE    COEFFICIENTS 

COr^KON/A/    GM(13?,36)*G(132i36)>GN(i32>36),S0(132)*Sl(132),S2(132) 

1  >A0( 132),A1(132)* A2(132  )» A3 (132 )» B0{36)>B1 (36)*B2(36) 

2  >B3(36),NX,NY*IX1,1X2,kSYM,FMACH,ALPHA,CA>SA,FMACH2 

3  >AL*UTIM*Ce>SB>NS,RG,IG,JG 
COMMON /B/  SV(132)»5r'(132)*CP(13  2) 

CO MM ON /C/  XP(260)fYP(2o0),Dl(26G),D2(2fc0)»D3(26G) 
COMMON/E/  CHORD* XM*CL^CD*CM 

COMMON/ L/  TCL(e01),TC0( 80I),TCM(b01),TCP(9,801),TCPS(9)*CLS,CDS 
1  *CMS*NITS*lJUHPfNSTEP*JSTEP>PERIOO*MHALF 

COMMON/WAKE/  NIT,WG(132)    •  ••' 
CL         =  0. 
CD         =  0. 
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12 


10 


11 


13 


1^ 


Ci-1         =0. 

N=IX2-1 

DO  12  I=IX1»N 

DX={XP{I  +  1)-XP(I)  )/CHJ(<D 

DY»  (YP( I+l  )-YP( I )  ) /CHORD 

XA=(.b*(XP(I+l)+XP(I))-XM)/CHQPD 

YA«.5*(YP(I  +  l)t-YP(I))/CHGRD 


CPA 

DCL 

DCD 

CL 

CD 

CM 

DCL 

CD 

CL 

IF(NIT.NE.O)  RETURN 

IF(NITS.EO.O)  GC  TO 

IJUMP=  2**MHALF 

CLS=  CL 

CD 

CM 

(  IX1  +  IX2  )*.i? 


+  CP  (I  )  ) 


.5*(CP( I+l) 

-CPAfQX 

C  PA*DY 

CL   iDCL 

CD   +DCD 

CM   +DCD*YA   -DCL*XA 

CL*COS( ALPHA)   -C D* S I N ( ALPh A ) 

CL*SIN(  ALPHA)   +CD*CCS(ALPHA) 

OCL 


11 


CDS  = 

CMS  = 

ISP  = 

1=  1 

TCPS(  I)=  CP( ISP) 

IFdSP.GE.  1X2  )  CD  TO  11 

1=  I  +  l 

ISP=  ISP+  IJUMP 

GO  TO  10 

NITS=  0 

JJSTtP=  JSTEP  +  1 

TCL(JJSTEP)=  CL-CLS 

TCD( JJSTEP)  =  CD  -CDS 

TCM(JJSTEP)  =  CM  -CMS 

IJUMP=  2**MHALF 

ISP=  (  IXI  +  1X2)*  .5 

1=  1 

TCP(  I^JJSTEP) =  CP(ISP) 

IF(  ISP. GE. 1X2)  GO  TO  1^ 

1=  I  +  l 

ISP=  ISP  +  IJUMP 

GO  TO  13 

RETURN 

END 


-rCPS  (I  ) 


SUBROUTINE  CPLOT 

PLOTS  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE 

COM/.ON/A/  GM(132,36),G(132>36)>GN(i32»36)^S0(132),Sl(132)>S2(132) 

1  >A0(132),Al(132),A2(132)>A3(132),B0(36),Bl(3o)>Bt;(3e) 

2  *33(36),NX,NY^IX1/IX2»KSYM,FMACH>ALPHA*CA*SA>FMACH2 

3  > AL*UTIM,CB»3b,NS* RG> IG* JG 
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12 


Z2 


falO 


COM 
CU'^ 
COM 
DI:^ 
DAT 
CPO 
CPO 
IF{ 

11  = 

12  = 

Wkl 
?  FOR 
1 
Du 
LIN 
DU 
K  = 
IF  ( 
IF( 
LIN 
WKI 
LIS 
kFT 
FDR 
END 


liDN/ 
MQN/ 
^ON/ 
FNil 
A 

IS 
=  0. 
ISTA 
IXl 
1X2 
IT 

TE  ( 

MAT( 

9 

12  I 

Ed) 

22  I 

20 

K.L 

K.G 

c  (K  ) 

Tfc  ( 

E(K) 

URN 

MAT( 


6/    5 V( 132)^  iM(132)>CP(132  ) 

C/  XP{260)*YP(2tjO)>Jl(2faO)>D2(2tO),D3(260) 

STADI/  PR*  IR» JR>  IbTADI 

UN   KGDE(2)^  LINE  (11!/) 

KODE/IH     >l.^+/ 
RESERVOR    PRESSURE    COEFFICIENT    WHERE    0=0    AND    F1T=0 

DI.oT.O)     CP0=((1.+.2*FMACH2)**3.5-1.)/(.7*FMACH2) 


=    5 
IWRIT^2) 

!)OHOPLnT  bF  CP  AT  tOuAL  INTERVALS  IN  THE  MAPPED  PLANE/ 
HO    X     »^H      CP   ) 
=  i»  115 

=  KuDE(  1) 
=11*12 

.♦ (CPO-CP(  1  )  )  +55.0 
T.l )   K=l 
T.115)  K=li5 

=  KDDE( 2) 
IWPIT^oiO) XP ( I ),CP ( I  )/LINE 

=  K  G  0  F  (  1  ) 

IH  ,2F9.^> 115Ai) 


SURRDUTINE   CHART 

GENERATES  MACH  ND  CHART 

CUMi^UN/A/    GM(i32>  36  )*G(  132*36  )»GN(132>36)*S0(132)*S1(132)>S2(1  32) 

1  >A0(132),Al(i32),A2(132)»A3(132)/B0(36)*Bl(36)*B2(36) 

2  /33(36)*NX*NY*IX1*IX2>KSYM*  FM ACH, AL PHA, C A, SA , FMACH2 

3  *  AL*UTIM*CB>S6>NS*RG,  IG*JG 
CQMMON/H/  DX, CY, DT*DXX,DYY,UTT,DXr*DXT,DrT,TSR 

CCMMQN/J/  RAC*PI*  ALS * AL T, AL TT > A MPL A » F RE QRA, F A  SAG  A/ FM ACHS * FMACHT 

1  *AMPLMiFRECRM>FASAGh>CETAS*CETAT>CETATT*AMPLC*FREQRC»FASAGC*CETA 

2  *FkEQR*IPSURE 

CDIiMCN/K/     10,  11,  12,  13,  Jl»  J2*  J3 

DIMENSION     IND(150) 

DDT=    l./DT  ■    ■  '  ■■  '        '  '      "  '  ■ 

AAO    =     1.  '"'■      ■"•■■  •  '       ■'•  '  '    ■  ■■  '    ■ 

I  WRIT  =6 

K  =    NY/32  '  '     ^      ■  '  '■  ^  ■   ) 

IF     (NY.GT.32*K)     K     =    K       +1 

WRITE    {  IWRIT, 2) 
2    FORMAT! l^HOMACH    NC    CHART)  -     '' 

11     DU    12     I «     10*  13  '.'.',.;.  ^ 

J  =  J1 

N=     0 
1^    N=    N+1 

X=     AC(I)  ^''    ' 
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12 


610 


Y=SO(  I)+BO( J ) 

Hh=  X*X  +  Y*Y 

DriH=  l./HH 

XT=  -. 5*Y*(ALT+CETAT)  +  UHH*(CA*X  +  SA*Y) 

YT=  .3*X*( ALT+CETAT  )  -DHh* ( C A*Y-S A*X ) 

H=SCRT{HH) 

DH«  l./H 

GI=«G(  I  +  1,J  )-G(I-l>J  ) 

GJ  =  G(  l» J  +  1)-G(I* J-1 ) 

GX         =  Aid  )*  GI   -SI  (I  )*B1  (J  )*  GJ 

GY         =   BKJ)*  GJ 

U=  &X*DH 

V=  GY*DH 

QQ=U*U+V*V 

CHAIN=  XT*GX  +  YT*bY 

FIT=  GM(  I, J  )  *DDT  +  CH4IN 

AA=  AAO  -.2*QQ-.4*HT 

AA=  ABS(AA) 

UR=  XT*H+  U 

VR=  YT*H  +  V 

QQR=  UR*UP  +  VR*VR 

QA=QQR/AA 

IND(N)=  1C0.*SQPT (QA) 

IF(UR.LT.O.)  IND(N)  =  -IND(K) 

J=  J+K 

IF  (J  .LE.  J2)  GO  TG  1^ 

WRITE  (I  WRIT, 610)  {IND(J),J'=1,N) 

RETURN 

FORMAT(  IX,  32U) 

END 


C 

c 
c 


1 

11 


SUBRGUT 
GENERAT 
AT  EQUA 
WITH  TH 
CDMMGN/ 

1 

2 

3 
CQMMGK/ 
COMMON/ 
COMMON/ 
COMMON/ 
COMMON/ 

1  ,AMPLM 

2  ,FREQR 
DIMENSI 
IF  (IPL 
CALL  PL 
CALL  PL 
CALL  SY 


INE  PSURE 

ES  PRE:>SURE  DlifRlBUTION  OVER  AIRFOIL 

L  INTERVALS  IN  THE  MAPPED  PLANE 

E  ASSOCIATED  SHUCKS 

A/  GM(i32,36),(.(132,36),GN(132,36),S0(132),Sl(132),S2(132) 
,A0(132)>A1(132),A2(132),A3(152),B0(36),B1(36),B2(36) 
,B3(36),NX,NY,IX1,IX2,KSYM,FMACH,ALPHA,CA,SA,FMACH2 
,AL,UTIM,CB,SB,NS,RG,IG,JG 

B/  SV(132  ),5M(132)*CP(132 ) 

C/  XP(2  6C),YP(2  6U),D1(26  0),D2(260),D3(26C) 

E/  CHaRD,XM,CL,CD,CM 

G/  TITLE(20), IPLOT 

J/     RAD,PI,ALS,ALT,ALTT,AMPLA,FREQRA,FASAGA,FMACHS>FMACHT 

»FPtQRM,FASAG(-;>CETAS,CETAT,CET4TT»AMPLCfFREQRC,FASAGC*CETA 

,  IPSURE 

ON     X(260),Y(  260)  ,R  (  IIjC  ) 
i»ll,101 


OT) 

OTSBL( 5000,23HI-CHUNG    CHANG 

0T(2.5,2.00,-3) 

MBuL(-2.0,-1.5-0,  .07,3,0.,-l  ) 


1C91C^W) 
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12 


l-V 


15 


22 


24 


25 


32 


34 


CA 
hU 
FG 
CA 
FA 
KA 
F  A 
EN 
FG 
1 
CA 
EN 
FO 
C  A 
EN 
FG 
CA 
XN 
XM 
DC 
XM 
Xi-l 

sc 

DG 
X( 
Y( 
N  = 
CA 
CP 
IN 
00 
AB 
IF 
CP 
IM 
CG 
C  A 
CA 
CP 
A  A 
CP 
IF 
ICA 
DG 
IF 
IF 
CA 
CD 
DG 
IF 
IF 
CA 
CG 
CA 


LL  SYMB 
CGDF( 43 
R  M  A  T  {  1 2 
LL  SYMB 
A=  FASA 
M=  FASA 
C=  FASA 
CODEC  57 
RMAT(5h 
F7. 
LL  SYMb 
CUCE (42 
RMAT( 5H 
LL  SYMB 
C3DE(42 
RMAT( 5H 
LL  SYi^b 
AX  =  XP(  I 
IN=XP(  I 
22  1=  I 
AX=AMAX 
IN 
ALE 

24  1  = 

I )=SCAL 
I)  =  SCAL 
IX2-IX 
LL  LINE 
NAX=  0. 
AX=  IXl 

25  1  = 
SCP=     AB 
( AeSCP. 
M  A  X  =  A  3  S 
AX=     I 
NTINUE 
LL    PLOT 
LL    AXIS 
C     IS    CR 
=     (1 .+. 
C= ( AA** 
(     CPC.G 
LL     SYMd 

32  1  = 
(CP(  I)  . 
(CPd  ). 
LL  SYMB 
NTINUE 

34  1  = 
( CP(  I)  . 
(CP(I  ). 
LL  SYMB 
NTINUE 
LL  SYMB 


GLCj 
>  12, 
A  4  ) 
0L(- 
GA*1 
G  M  *  1 
GC*1 
,  14, 

1  IMF 

2  ) 

GL(- 
,15, 
AL 

UL  (- 
,16, 
CL 

QLC- 
XI) 
XI) 
XI,  I 
1  (XP 

=  AN 
=  5 
1X1, 
£*(X 
t*YP 
1+1 
(X(I 


1X1, 
S(CP 
LE.C 
CP 


(C, 
(-1. 
ITIC 

2  *  F  M 
3.5- 
E.-l 
GL(- 
IXl, 
GT.l 
LT.- 
GL  (X 

IMAX 
GT.l 
LT.- 
QL(X 


.5,-1. 50, .07,3,0. ,-1  ) 
R)  TITLE 

.5, -.75,  .14,R,0. ,48) 

80. /PI 

8  0  .  /  P  1 

80. /PI 

R)  UTIM,FAA,FAM,FAC 

=  ,F7.2,3X,jHALFA=,F7.2,3X,5HM  FA'=,F7.2,3X^5HCEFA», 

.5,-1.0,.14,R,0.,5  7) 

R)  AL,FMACH,CFTA 

=  ,F7.3,3x,5H,-l    =,F7.3,3X,5HCETA-,F7.3) 

,5,-1.  25,  .14,K,G.,42) 

R)  CL,CD,CM 

=  F7.4,3X,5HCD   =  ,  F 7  .  4, 3X , 5HCM   =,F7.4) 

.5,-1.50, .14,R,0.,42) 


X2 

( 1), XMAX) 

INI (XP ( I),XMIN) 

./(XMAX   -XMIN) 

1X2 

P( I )-XMlW) 

(  I) 

Xl),Y(IXl),N,l,G,i,0.,l.,C.,l.) 


1X2 
(  I)  ) 

PMAX )  GG  TG  25 


4.25,-3) 

,-4.,2HCP,2,8.,90.,l,6,-.4,0) 

AL  PRESSURE  CGEFFICIENT 

ACf;2)/1.2 

!.)/(. 7+FMACH2 ) 

.6.AND.CPC.LE.1.6) 

l.,-2.50»CPC, .4, 15,0.,-1) 

IMAX 

.6)  GO  TO  32 

l.t)  GO  TG  32 

(I),-2.5  0*CP(I),.07,3,45.,-1) 

,1X2 

.fc)  GO  TO  34  .  .  ..  , 

l.t)    GO  TJ  34 

(I),-2.5  0*CP(I),.07,3,0.,-1) 


0L(-2.0>-5.75,.07,3/0.,-l) 
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101 


CALL  SYMBDL(6.l5»-b.7i<*.07,i,0./-l) 

CALL  PLGT  (-2.5>-6.2tJ/-3) 

CALL  FRAME(l) 

RETURN 

CALL  PLaT(0.>C.^999) 

RETURN 

END 


C 

C 


12 


14 


15 


le 


SUBK 
GENE 
AT  E 
COMM 

1 

2 

3 
CQ.-^M 
CGMM 
CGMM 
CQfiM 
COMM 
CD^^M 

1  ,AM 

2  ,FP 
DIME 
11  =  1 
1^  =  1 
IF  ( 
CALL 
CALL 
CALL 
CALL 
ENCO 
EDRM 
CALL 
I-  AA  = 
EAM  = 
FAC» 
ENCO 
FORM 

1 
CALL 
ENCO 
FORM 
CALL 
ENCO 
FORM 
CALL 
XMAX 
XMIN 
DO  2 
XMAX 


XP  ( 
CHG 
XR> 
TIT 
RAD 


DUTINE  LO 

PATES  THE 

OUAL  INTE 

ON/A/  GM( 

^A0( 

*  B3( 

,  AL» 

ON/e/  sv( 

UN/C/ 
UN/E/ 
DN/F/ 
Ci  N  /  G  / 
ON/J  / 
PL,1,FKEJK 
EGR, IPSUR 
NSIQN  X{2 
XI 
X2 

IPLGT)  1, 

PLOTSBL( 

P  L  G  T  (  2  .  5 

SYM3DL(- 

SYMBuL  (o 

L;E(  46>  12, 

AT( 12A4) 

SYMBGL(- 

FASAGA*! 

FASAGM*1 

FASAGC*! 

DE(37»14, 

AT( SHTIME 

F7.2) 

SYMeQL(- 

DE(42, 15, 

AT(5HAL 

SYMBOL (- 
0E(42,lb, 
AT(5HCL 

SYMBDL(- 
=  XP(  ID 
=XP(I1) 
2     1=I1»I2 
=AMAX1 (XP 


RO 

LOADING    DISTRIBUTION    OVER    AIRFOIL 
RVALS     IN     THE     MAPPED    PLANE 

13?,36),G(132,36),GN(132,36),S0(132),S1(132),S2(132) 
132),A1{132),A2(132),A3(132),B0(36),B1(36),B2(36) 
36)>NX,KY,IX1,IX2,KSYM,FMACH,ALPHA,CA,SA,FMACH2 
UTIM,CB»S3,NS,kG,  IG> JG 
13?),SM(132),CP(132) 

260),YP{2  6O),Dl(260),U2(2t)0),D3(260) 
RD,xr,  CL,CD>Cr^ 
YR,KS, XS (500) ,YS  (  500) 
LE (20), IPLGT 

,PI,ALS,ALT,ALTT>AhPLA,FKEQ-^A,FASAGA,FMAChS,FMACHT 
M,FASAGM,CETAS,CEIAT,LETATT,AMPLC,FREOftC,FAbAGC>CETA 
E 
bo)  ,  Y(2t.O),K(150),DCP(  132) 


11,101 

5C00,23HI-CHUNG    CHANG  1C9104W) 

,2.00,-3) 

2.C,-1 .50, .07,3,0.,-! ) 

.5,-1.5C,.u7,3,C.,-l) 

R)     TITLE 

.5,-.75,.14,R,0.,'fH) 

80. /PI 

dO. /PI 

60./ PI 

R)     GTIM,FAA,FAM,FAC 

=,F7.2,3X,5HALFA=,F7.2,3X,5HM     FA=,F7.2,3X,5HCtFA», 

.5,-1.0, .14,R,0.,57) 

R)     AL,Fr^,ACH,CETA 

=',F7.3,3X,5HM  =,F7.3,3X,5HCETA«,F7.3) 

.5,-1.25, .14,R,C.,42) 

R)    CL, CD, CM 

=  F7.'t,  3X,  5HC0       =,F7.'V,3X,5HCM       =,F7.4) 

.5,-1.50,. 14,R,C  .,42  ) 


(  I  ), XMAX) 
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22  XKIN        =A,^Iin  (  Xf^  (  1  JiXMIK  ) 

SCALE       =  5./(XMAX   -Xi1IN> 

DG  2^  1=  11*  12 

X(  I  )=SCALE*(XP(  I  )-XMI:M) 
24     Yd  )  =  i>CALE*YP  (I  ) 

IN3SF=     .3*{  Kl    +     IX?) 

Du  26  1=  INUSE^IXZ 

>,=  ro  +  4  -  I 

DCP(I)   =   CP(M)-CP(I) 
2&  CGNTINUt 

N  =  I?   -II   +1 

CALL  LINt  (X(Ii)>Y{Il)/;J,i>C^l^C.^i.*0.^1.) 

CALL  PL0T(0.,4.2!5^-3) 

CALL  AXIS(-l.*-4.>3HDC?.3^8.,9C.,-l.fc,.4,0 

DO  32  I  =  INJSE,'I2 

It- (  JCP  (  I  )  .GT.  1.6)  GO  T[j  32 

IF (OCP( I) .LT.-1.6)  GC  TO  32 

CALL  SYM80L(X{I),  2 . 5*DC P ( I ) » .0 7 , 3^ 0 . , -1 ) 
32  CONTINUE 

CALL  jYMbGL(-2.C,-:/.7:>*  .O7*3»0, 

CALL  SYMBOL  (t  .5,-5.75/ .07,  3, O.J 

CALL  PLQT (-2. 5,-6.25,-3 ) 

CALL  FkAME (1) 

RETUkK 
101  CALL  PLOT (0. ,0., 999  ) 

RETUKN 

END 


,-1) 
-1  ) 


SOBRCL'TINE  SONIC 
C      GENERATES  SONIC  LINE  OVER  AlkFCIL 
C       RENAME  COMMON/A/ 
C      THE  POSITION  OF  GN  IS  OVERLAPPED  BY  SHOCK 

CLMMON/A/  GM( 132, 36),G(132,36),SH0CK(i32,36) 
*  ,S0( 132), Sl( 132), S2( 132) 

1  ,A0(132),Al(132),A2(132),A3(132),30(36),Bl(36),fl2{36) 

2  ,B3(  3e  )  ,NX,NY,  IXl,  IX2,KSYM,F,-iACH,ALPHA,CA,SA,FMACH2 

3  ,  A  L  ,  U  T  I  M  ,  C  B  ,  S  3  ,  N  S  ,  ??  G  ,  I  G  ,  J  G  ^  .   .,,..... 
CuMMGN/B/  SV(  i32),SM(132),CP(132)  .•  ...; 
COMMGN/C/  XP( 260) ,YP(260), 01(260), 02(260), 03(260)        '-  -k 
COMMON/D/  SLOPT, TRAIL, SCAL    .  y  »..,„.,(.,  i ..  /  -  ,   ., 
CuMMGN/E/  CHORD,XM,CL,CD,CM 
COMi"lGN/F/  XR,  YR,KS,XS(500),YS  (500) 

COMMLN/G/  TITLE(20),  IPLOT  ,..,.,  t  ^ .'  ,  ..  , 

COM.-iCN/H/  OX,  DY,DT,DXX,DYY,DTT,DXY,DXT,DYT,TSR 

COMMON/J/  RAD,PI,ALS,ALT,ALTT,AMPLA,FREuRA,FASAGA,FMACHS,FMACHT 

1  ,AMPLM,FREQkM,FASAGM,CETAS,CETAT,CETATT,AMPLC,FREQRC,FASAGC,CETA 

2  ,FREQR,IPSURE 

CCMMON/K/  10,  II,  12,  13, J1,J2,  J3 

DIMENSION  XA(2bO),YA(2oO),k(150)  , 

DIMENSION  C{260) 

IF  (IPLOT)  1,11,101  ^ 


\t.)^-  ->:  ■* 
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1 

11 


12 


14 


lb 


16 


22: 


24 


CALL 
CALL 
CALL 
CALL 
ENCGD 
FORMA 
CALL 
FAA  = 
FAM  = 
FAC  = 
ENCDD 
FORMA 
L 

CALL 
ENCnD 
FLiRMA 
CALL 
FNCOC 
FORMA 
CALL 
AIRFC 
X  f*i  A  X  = 
XMIN  = 
DO  ZZ 
XMAX  = 
AhIN 
SCALE 
DU  24 
XA(  1) 
YA(I) 
CALL 
N=  IX 
CALL 
AAO  = 
DDT  = 
DO  2 
DO  2 
X=  AC 
Y=SO{ 
riH=  X 
DHH  = 
XT=  - 
YT=  . 
H=SJk 
0H=  1 
GI^GC 
GJ  =  G( 
GX 
GY 

U=  oX 
V=  GY 
UU  =  U* 
VV  =  V  * 


PLOTSBL  (  i?000.  23HI-CriUNG  CHANG    lO^lO^W) 

PLOT( 2.5,2.00,-3) 

SYMB0L(-2.G,-1.50>.07,3>0.,-1) 

SYMBOL (6. 5,-1. 50,. 0  7, 3, 0.,-l) 

E(48,12,R)  TITLE 

T(12A4) 

SYMB0L(-.5,-.75,.i4,K,0.,46) 

FASAGA*180./PI 

FASAGM*180./ PI 

FASAGC*180. /PI 

E(57,14,R)     UT1K,FAA,FAM, FAC 

T(5HTIMe=,F7.2,3X,5HALFA=,F7.2,3X,5HN     FA=,F7.2,3X,5HCEFA=, 

F7.2) 
SYMB0L(-.5,-1.0,.l4,R,0.»5  7) 
E(42,15,R)     AL,FMACH,CETA 

T(5HAL        =,F7.3,3X,5HM  = , F 7 . 3 , 3 X, 5HC E T A  =  , F 7 . 3  ) 

SYMB0L(-.5,-1.25,.14,F,0.,42) 
E(42,lb,R)     CL,CO,CM 

T(5HCL   =F7. 4, 3X,5HC0   = , F 7. 4 , 3 X, 5HCM   =,F7.4) 
SYMBOL (-.5,-1. 5C, . 14, P,Q., 42) 
IL 
XP(  1X1  ) 
XPCIXl) 
1=  1X1,1X2 
AMAXl (XP  (  I  ), XMAX) 

=AMIN1(XP( I ) ,XMIN) 
=  5. /(XMAX  -XMIN) 
I =  1X1,1X2 

F  * (aP( I )-XMIN) 
*YP(  I  ) 
,4. ,-3) 


PL 
2- 
LI 
1. 
1. 
1  = 
J  = 
(I 
I) 
*X 
1. 
.5 
5* 
T( 
./ 
1  + 
I, 


*D 
*D 
U 
V 
+  V 


=   IXl; 

SCALI 
SCALE' 
GT(U. 
IXl+1 
f4fc(XA(IXl),YA(Ixl),N,l,0,l,0.,I.,0.,l.) 


/UT 

10 
Jl 

) 

+  B0 
+ 

/HH 

*Y* 

X*  ( 

HH) 
H 

1,J 
J+1 


,  13 
,  J3 

(J) 
Y*Y 

(ALT+CETAT)  +  DHH*(CA*X  +  SA+Y) 
ALT  +  CETAT)  -IJHH*  (  C  A*  Y-SA*  X  ) 


)-G( I-l, J ) 

)-G( I, J-1) 

=  Al( 1 )*  GI   -Sl( I)*f 1( J )*  GJ 

=   Bl( J  )*  GJ 
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IQ 


20 


le 


21 


52 
70 


101 


'Jk=     XT*H  + 
y/  R  =     Y  T  *  H 
UU-<=    UR*U 
VVR=    VR*V 
&UR=    (JUR 
CHAIN=    XT 
FI r=     GM( 1 
AA=  A/>0-.2 
AA  =  Ai'^AXl  ( 
bhJCK(I>  J 
CCjNTINUE 
LLCATES     T 
•<>S=    0 
L^J     17     J  = 

w( J)     =    Sh 

DO   la   1= 

Du  It  K= 
J=  J1+J2- 
wv,=  SHOCK 
I  F  (  a  Q  .  N  E  . 

K  J  =   K  5  +1 

XS(KS)  = 
Y  S  ( .<  S  )  = 
GLi  TG  16 
IF ( J&«G( J 
f.T=  AbS(wC 
K3=  KS  +1 
Xb  (KS  )   * 
Y1,(KS  )=  S 
IF (CU*Q ( J 
?T=ABS(  Qi. 
KS=  KS  + 
XS ( KS  )=  A 
YS(KS) =S0 
u!(  J  )=  QQ 
Ou  21  1= 
XX=  .5*SC 
YS(I  )  =  J 
X  S  (  I )  =  X 
OG  :;?  1  = 
XX=  SCALE 
YY=  SCALE 
IF  {XX.  LT. 
IF(XX.GT. 
IF( AbS( YY 
CALL  SYMB 
CONTINUE 
CALL  SYMB 
CALL  SYMB 
CALL  PLOT 
CALL  FRAM 
RETURN 
CALL  PLOT 
RETURN 
END 


U 
+     V 
k 
R 
+    VV 

*r,x 

.J  ) 
*UQ- 
AA»  . 
)     = 


+     YT*GY 

*DCT     +    CHAIN 

.^*FIT 

0001  ) 

S  0  K  T  {  Q  0  P  /  A  A  )     -    1  . 


He     SONIC     POINTS 


J1»J 

LCK( 

II,  I 

J1>J 

K 

(  If  J 

0.  ) 

A0(  1 
S0(  I 

+  1  )  . 
/  { u;. 

A0(  I 
0(1) 
)  .GE 
/(Ow 

1 

0(1) 
(1  )  + 

i>KS 

AL*( 

CAL* 

X 

l^KS 

*(XS 

*YS( 

-2.0 

fc.  5) 

)  .GT 

CL(  X 

0L(- 
0L(6 
(-2. 
b  (1) 


3 

10^  J  ) 
2 

) 

GO    TO     19 


) 

)     +     B0(J) 

GE.O  .  )     GJ    TC    20 
-C(  J+1)  )  ) 


) 

+  BC(J)     +RT*(BG(  J  +  1)-cjO(  J  )  ) 

.0.  )    GO    ro    Id 

-0  (  J  )  )  ) 

+RT*( A0( I-1)-A0( I ) ) 
riO (J  )  +  RT*(SO( I-l )-SO(I  )  ) 


XS ( I  )**2-YS ( I )**2  )     +    XR 
XS(I)*YS(I)     +    YR       ,,- 

{  I  )-XMlN) -.   ^,   ^      . 

I) 

)     GO     Tu     52  ,         .     0-     • . 

GO     TO     52 
.4.5)    GO    TO    52 
X, YY, ,07, i,0.  ,-l) 

2.0,-5.5    , .G7,3,0.,-l) 
.5,-5.5    ,.07,3,0.,-!) 
5,-6.0,-3) 


{0.,0.,999) 
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1 

11 


12 


lA 


15 


lb 


SUBRCUTINc  TkACE 

GENERATES  UNSTEADY  TRACES  UP  AIPFGIL 

COMMON/A/    Gr(132,3D)»G(li^»36)»GN(i32,36),S0(132)>Sl(132)»S2(132) 

1  ,AO(132)^Al(132)>A2(132)^A3(132),BO{36),Bl(36)»B2(3t) 

2  ,B3( 3b )*NX*NY*IX1*IX2*KSYM»FMACH, ALPHA, CA>SA*FMACH2 

3  ,  AL,  UTI^1KB,SB,iMS,kG*IG,  JG 
REMAME     CUi^MUN/C/ 

CUMMCN/C/     XP(2  60),YP(2o0),Dl(2b0),Y(2  6C),R(2  60) 

CGMMuN/G/     TITLE(20)  WPLCT 

CDMMUN/H/     0X,DY*DT,DXX>DYY>DTT>DXY,DXT*DYT,TSR 

CGMMCN/J/     RAu,PI>ALS>An,ALTT,Ai'«'PLA,FREOKA,FASAGA,FMACHS,FMACHT 

1  ,AMPLM>FREaRM,FASAGM,CETA5,CETAT,CETATT>AhPLC,FREQRC>FASAGC,CETA 

2  iFREQR,  IPSURE 

CQMMCN/L/     TCL(8cl),TCD(801),TCt'',  (SDD^TCPCy^aOD^TCPSCQj^CLS^CDS 
1  ,CMS»NITj,lJU-^P*NSTEP*JSTEP,PERIOD>MHALF 

DIMENSION     X(H01),AuA(601)/FS(80i),FA(8Cl) 
IF     (iPLUT)     1,11*101 

CALL  PLGTiBLC 5CCC,23hi-CHUNG  ChAMG    IC^IO^W) 
CALL  PLOT( 2.5,2,00,-3) 
CALL  5YMBDL(-2.C,-;..5C,.0  7,3,0.,-1) 
CALL  SYMBuL(6.i',-1.50,.07,3,C.,-l) 
TITLE 

cNCLlLE(^b,12,  R)  TITLE 
FuKMAT( 12A^) 

CALL  SY-1BUL(-.5,-.5C,.i't,R,0.,'^o) 
LNCGUt (SfiflH,  R  ) 

rURMAKACHUNSTEADY  TPACES  OF  AIRFOIL  IN  SINUSOIDAL, 
1        18H  RIGID  BODY  MOTiON) 
CALL  SYMB(JL{-.5,-.75,  .1^,R  ,0.,'j8) 
ACAS=  ALS*RAD 

tNCQCE(57,15,R)AGAS,A'^PLA,FRtORA 
FGRriAT(  leHMEAK  ATTACK  ANGL  t  =,  F  5  .  2,  5  X  ,  ^HAMP  =  ,  F5  .  2  ,  5  X  , 

lOhFREO  ^ATE=,F5.2) 
CALL  SYN6UL(-.5,-1.0C,  .1^,K,0., 57) 
ENCGD£(57,i6,P)  FMAChS,ArlPL.'^,FPEQRM 
FORMATdBHKt  AN  FLIGHT  S  Pb  t  D  =  ,  F5  .  2,  5  X,  <tH  AM  P=  ,  F  5  .  2  ,  5X, 

lOHFRFQ  RATE=,F3.2) 

1^, K, 0., 57  ) 


1 


CALL  SYii80L(-.5,-l  .25, 

FANGLE=  CETAS*RAD 

ENCODE (37, 17, k)  F ANGLE, AM  PLC, FREQRC 
17  FURMAT( IBHMEAM  FLIGHT  ANGLE =, F 5 . 2, 5 X , ^HAMP « , F 5 . 2, 5X , 
1       lOHFkEu  RATE=.F5.2) 

CALL  SYMBL:L(-.5,-1.50,.l'^,R,0.,57) 

AIRFOIL 

II  =  IXl 

12=  1X2 

XMAX  =  XP ( II  ) 

XhIN=XP(  ID 

DQ  52  1=  II, 12 

XMAX= AMAXl  (XP  (  I  )  ,XMAX  ) 
52  XMIN       = AMINl ( XP( I ),XMIN) 

SCALE      =  3./(XMAX   -XMIN) 

DO  5^  1=11,12 

X(  I)=SCALE*(XP(1  )-XMIN) 
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14 


^C 


J  'J 
5b 


b7 

5P 


59 
bC 


y  (I )  = 

In 

CALL 

CALL 

u  K  A  W  b 

IJ  jrP 

I5P  = 

IPQIN 

CALL 

CALL 

hNCub 

FuR;1A 

CALL 

I  K  (  I  P 

I  F(IS 

IPGIN 

I^P  = 

Gu  TU 

UNSTE 

XSCAL 

IF  {L^' 

ASCAL 

Gu  r  L 

ASCAL 

IF{  A^ 

FSCAL 

Gu  TC 

F5C  AL 

I  F  (  A 1'^' 

CSCaL 

&L]  TC 

CSCAL 

JS=  J 

DL  2 

T  I  M  t  = 

X(  I  )  = 

ALA(  I 

FS(  I  ) 

f-  A(  I  ) 

CGNri 

CALL 

CALL 

EnClJU 

F  U  P  M  A 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

DG  3 

CALL 

CQ.Ul 


SCALF->YP(  I  ) 

=12   -II   +1 
^  .  »  -  .  ^  :;  »  -  3  ) 

SUft  SENSuRS  UN  THE  AIRFuIL 

^HALF 

1X2  )*.5 


PLO 
LIN 

OR 

=  2 
(  IX 
T- 

SYM 
SYM 
E  (1 

r(i 

SYM 
GIN 

P.o 

T  = 

ISP 
19 
AOY 
=  ( 
PLA 
=  .2 
56 
=  0 
PLrl 

53 
=  0 
PLC 

"    • 

60 
-  0 
STE 
1=1 
DT 
XS 
)  = 
=  F 
=  CS 
NUE 
PLQ 
AX  I 
t  (  i 
T(l 
SYM 
PLD 
PLU 
PLC 
SYM 
SYM 
PLQ 
I« 
PLD 
KUE 


T  ( 

E  ( 
£S 
*« 

1  + 
1 

dC 
QG 
,  1 
1  ) 
BC 
T. 
i  • 
IP 
+ 


L (X  ( ISP )*  Y (  liP)^  .07,3,0.^-1) 
L(X(ISP),Y(ISP),.C7,3,45.,-1) 

a  ,  K  )  I  P  D  1  N  T 

L( X(  I SP ) ,  .1  +  Y (I SP  ),  .G7,k,0.,l ) 
Gc . 9)  GJ  Tu  20 
1X2)  GG  TC  2u 

QINT  +  1 
IJ  U  (^  P 


PkEiSukE     TRACES    ON    Tht     PRESSURE 
e.*FR£CR )/ (2.*PI*PtKlCD) 
. EG.O. )    GG    TG    5  5 
/AMPL A 


SENSORS 


.  E  0  . 0  .  ) 
2  /  A  M  p  L  >■- 


.EQ.C. ) 
2/ AT: PLC 


G    TD    5  7 


-C    TQ    59 


P     +    1 
,  JS 
*NSTEP*( I-l ) 

C  A  L  *  T  I M  E 

ASCAL*Ah,  PLA*SKaTIME*FKEQRA) 

SCAL*AK,PLM*SiN(FKECRr*TIME) 

CAL*AMPLC*SlN(FRf.iRC*TIME) 


T{-2.  5,  .75,-3) 

S  (  0  .  ,  0 .  ,  1 1 H 

100,R) 

IHPHASEANGLL) 

BUL(0.,-.5,.14,R,0.,11) 

T(     0. ,0., 3)  ,     ,    . 

T  (  0  .  »  b  .  ,  2  ) 

T (0., . 5,-3) 

aOL  (CO.,  .07,15,0.,-1) 

B0L(-.75,0.,.1^,6HAANGLE,0.,6) 

T(     0.,C., 3) 

1,  JS 

T(X( I  ), AGA(  I), 2) 


,-ll,6.,0.,0.,.b0*PERI0D,0) 


.'  t 


I  (  ■  1 


J  f  r 
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CALL  PLOT(0.>  .5,-3) 

CALL  SYMBQL{O...C.,  .G7,15>0.,-1) 

CALL  SYMdLLC-.  75,0.,.  l-t^bH  FAN  GLE>0.»6) 

CALL  PLQT(  0.,0.,3) 

DU  ^  1=  1, JS 

CALL  PLOT (X(I  ),FA (I ),Z) 

CONTINUE 

CALL  PLOT( 0., .5,-3) 

SYMBOL  (CO.,.  07,15,0.,-!) 
SYMBOL(-.7i;,0.,.i^,6HFSPEED,0.,fc) 


71 


7o 
7^ 


72 


CALL 
CALL 
CALL 
DO  5 
CALL 


73 


CALL 
CALL 
CALL 
00  6 
CALL 


5,-3) 

,0  .  ,  . 
26,  C. 


07,15,0.,-!) 

,  .  I't  ,2MC^. »  0.  ,  2  ) 


PLOT (  0.  ,0.,3) 

1=  !,JS 

PLOT(y ( I ),FS ( I),2) 
CONTINUE 

IF  (KSYM.GT  .0.  ANIJ,  ALS.n- 
TCMMAX^     0. 
DO     7!     1=     !,Jb 
ABSTCM=     ABSCTC^d)) 
1F( ABSTCM.LE.TCMMAX)     oG     TO    71 
TCi1MAX=     ABSTC^ 
CONTINUE 

TCMCAL=     .2/TCMMAX 
GO    TO    7^ 
TCMCAL=0. 
CALL     PLQT(0., 

SYMBOL (0 

SYMBOL (- 

PLOT(     O.,0.,3) 

1=     1,JS 

PLGT(X( I ), TCMCAL*TCM( I), 2) 
CuniNUE 
TCDMAX=  0. 
DC  72  1=  1,JS 
ABSTCD=  AdS ( T C 0 (  1  )  ) 
IF( ABSTCD.LE.TCLMAX)  GO  TO  72 
TCDMAX=  ABSTCO 
CONTINUE 

TCDCAL=  .2/TCDfAX 
CALL  PLOT(0.,  .5^-3) 
CALL  SYMBOL  (CO.,  .07,15,0.,-!  ) 
CALL  SYMB0L(-.2t,0.,  .14,2hCL',  0.  ,2) 
CALL  PLaT(  0.,C.,3) 
DO  7  1=  !,JS 

CALL  PLOT(X(  I  ),TCDCAL*TCD(  I)/2) 
CONTINUE 

IF(KSYM.GT .0. ANO.ALS 
TCLMAX=  0. 
DO  73  1=  1,JS 
ABSTCL=  ABS(TCL  (1  )  ) 
IF(A«STCL.LE.TCLf.AX) 
TCLMAX=  ABSTCL 
CONTINUE 

TCLCAL=  .2/TCLMAX 
GO  TO  75 


AND.AMPLA.tC.O..AND.A^iPLC.tQ.C.)  GuT076 


EJ.G..ANb.AMPLA.EC,O..AND.AMPLC.EC.O.)  GCT077 


GO  TO  73 
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77 
Id 


3 'J 


?.•=> 


'^l 


A2 


^3 


TCLC 
CALL 
CALL 
CALL 
CALL 
Du  a 
CALL 
CUNT 
TCP'1 
DL  3 
OU  3 
ABST 
IF(  A 
TCPM 
Cl.i4T 
TCPC 
CALL 
CALL 
CALL 
CALL 
\)u  9 
CALL 
CJNT 
CALL 
CALL 
CALL 
CALL 
DC  2 
CALL 
CUNT 
CALL 
CALL 
CALL 
CALL 
Du  4 
CALL 
CONT 
CALL 
CALL 
CALL 
CALL 
DU  4 
CALL 
CUNT 
CALL 
CALL 
CALL 
CALL 
00  ^ 
CALL 
CUNT 
CALL 
CALL 
CALL 


AL  =  0 
PLU 
SY;i 
SYM 
PLO 
1  = 
PLG 

INUE 

AX  = 

C  K  = 

C  1  = 

CP  = 

bSTC 

AX  = 

lUUE 

AL  = 
PLG 
SY»1 
SYM 
PLC 
I* 
PLU 

I'VJE 
PLU 
SYM 
SYM 
PLO 

9  1  = 
PLO 

INUE 
PLU 
SYM 
SYfl 
PLO 

1  1  = 
PLO 

INUE 
PLU 
SYM 
SYM 
PLO 

2  1  = 
PLO 

INUE 
PLO 
SYM 
SYM 
PLO 

3  1  = 
PLO 

INUE 
PLG 
SYM 
SYM 


,  .01 , i  5>0.> -1 ) 

G.  »  .  1-^><::HCL>0.  f  2) 

3) 


T(0.^ .5»-3) 

b  C  L  (  0  .  >  0  . 

e  U  L  (  -  .  2  6  f 

T (     u. ,0., 

U  JS 

T(X(I  )^1CLCAL*TCl (1)^2) 


0. 

i»  JS 
AL>3  (TCP  (K 
P.  LF.TC  PT'i 
A3STCP 


.2/TCP>IAX 
T  (0.,  .b,- 
eUL (G.,0. 
qGL (-. 28, 
T  (  0  .  >  C  .  , 
1 ,.  J  5 
T ( <(I  )  ,TC 


>  I  )  ) 

AX)     Gu    TO    3C 


3) 

t  .C7, 13,0 . ,-1) 

0. , .  1^, 2HP1,0.,2) 

3) 

PCAL*TCP ( 1,1 ),2 ) 


3) 


T  (  G  . »  .  'j ,  - 
rt  U  L  (  0  .  ,  0  . 
iiuL  (  -  .  2t) , 
T  (     0 .  ,  C  .  , 

1»  JS 
T(X( I  ),TCPCAL*TCP( 2,1 ), 2) 


,  .07,15,0.,-! ) 
0.,  .1^,2HP2,0., 2) 
3) 


•3) 


T  (  G  .,  .bf- 
B  U  L  (  0  ,  ,  C  . 
DuL(-.2e, 
T (     0. ,0. , 

1,  JS 
T (X( i  ),TCPCAL*TCP( 3,1  ), 2) 


, .07,15,0.,-!) 
0.,  .i^,2HP3,0. ,2) 
3) 


■3) 


T(0., .5,- 
6  0  L  (  0  .  ,  C  . 
B0L(-.2fc, 
T (     0. ,0., 

!,  JS 
T{X(I  ) ,TCPCAL*TCP( 4,1 ) ,2) 


, .07,15,0.,-!  ) 
0,,  .l'^,2HP4,C.,2) 
3) 


3  )  .  :    ,  . 

,.07,15,0.,-!) 
0., .14,2HP5,0.,2) 
3) 


T  (  0  .  ,  .  5, - 
30L(0.,0. 
BUL(-.2b, 
T(     0. ,0., 

!,J5 
T (X( I ),  TCPCAL*TCP( 5,1 ),2) 


T(0.,  .5,- 

BC.L(G.,U. 
BDL(-.26, 


3) 

,  .07,15,0. ,-!) 

0.,  .  i't,2HP6,0.,2  ) 


■;  t 
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CALL  PLOT(  0.^0., 3) 
DO  ^^  1=  1^ JS 

call  plm(x(i)»tcpcal*tcp(6»i)»2) 
^a'  cuntinuf 

call  plot (0.> .5,-3) 

CALL  SY(13L1L(0.,0.  ,.  07,15  »0.,-l) 

CALL  SYM6LL(-.26,C.,.i^,2HP7,0.,2) 

CALL  PLQT(  0.»C.,3) 

UU  <t5  1=  1,JS 

CALL  PLGT(X(I),TCPCAL'!'TCP(7,I),^) 

45  CONTINUE 

CALL  PL3T(C., . 5,-3) 

CALL  SYMBQL(0.,0.*.0  7,15,0.,-i> 

CALL  SYMBGL(-.2b,0.,.14,2HPb,O.,2) 

CALL  PLOTl  0.,C.,3) 

DO  46  1==  1,JS 

CALL  PLGT(X(I),TCPCAL*TCP(8,I),2) 

46  CONTINUE 

CALL  PLOTCO.,  .5,-3) 
CALL  SYMBOL(0  .,0.,.G?,i5,0.,-l) 
CALL  SYMB0L(--.2b,0.,  .14,2HPv,u.,2) 
CALL  PLOT(  CO. ,3) 
DO  47  1=  1,JS 

CALL  PLar(X(l),TCPCAL*TCP(9,l  ),2) 
4  7  CONT INUE 

CALL  PLOT (  .5»-8.00/-3) 

SYMROL(-2.0,-i.50,.07,3,0.,-l) 

SYMBOL (6.5,-1.50, 

PLOT (-2.5.-2.0,-3) 

FRAME( 1) 
HETUPN 
101  CALL  PL0T(0.,0.,<*99) 
PETURN 
END 


CALL 
CALL 
CALL 
CALL 


07, 3,0.,-l) 


SOBKGUTINE     GRID 
C  PLOTS     THE     KtSH    bYSTE^. 

C  RENAME     COMflON/A/ 

C      THE  POSITION  OF  GM  AND  GN  ARE  OVERLAPPEO  bY  XhESH  AND  YMESH. 
C      THE  ROUTINE  SHOULD  BE  lALLED  AT  THE  kIGHT  END  OF  THE  PROGRAM 
COMMON /A/  XM£SH(132,36),b(132,3fc),YMESh(13::,36) 

*  ,S0( 132),S1 (132 ),S2(132) 

1  ,A0(13^),Al(i32),A2(132),A3(132),B0(36),61(36),62(3fc) 

Z  ,B3(3fc),NX,NY,IXl,IX2»KSYM,FMACH,ALPHA,CA,SA,FMACH2 

3  , AL, UT1M,CB,33,NS,PC, lb, JG 

commqn/d/  sljpt, trail, SCAL 
common/f/  xr,yr,ks,xs(300),ys(50c) 

COMMON/G/  TITLE (20), IPLOT 
COMMON/K/  10,  II,  12, 13, J1,J2» J3 
IF  (IPLOT)  1,11,101 
1  CALL  PLOTSBL( 5000,?3HI-CHUNG  CHANG    1C9104W) 


181 


11 


12 


1^ 


13 


33 
32 


-^3 
42 


101 


ALL 
ALL 
ALL 
NCOD 

(JRMA 

ALL 

NCJD 

0>^MA 

ALL 

ALL 

tSH 

0=  X 

0=  Y 

J  13 

G  13 

n  E  S  h 

MESH 

KAWS 

MAX» 

MIN  = 

G  22 

MAX  = 

M  1  N  = 

CALE 

G  3  2 

P 

C  32 

P=  S 

P=  S 

F(  XP 

ALL 

P 

U  TIJ 

P 

CNTI 

RAWS 

L  42 

P 

U  M 

P'  S 

P=  S 

F(  XP 

ALL 

P 

0  rc 
p 

ONTl 

ALL 

ALL 

ALL 

ALL 

ALL 

tTUR 

ALL 


PLDT(  2.  i;^2.0,-3) 

SYMBOL  (-2.  U>-I.!)0i.07>3^0.»-1) 

SY-^BLL{6.b»-l.t'0,.0  7f3>0.»-l) 

E(80>12,R)  riTLt 

T(20A4) 

SYMBuL( 0.^-. b^ .14,R>0., fcO) 

E(35/14,^^)  NX, NY 

T(24hNEAR  f-IELC  GRID  SYSTEM 

SYi'-'eDL(0.,-.75,.l4,R,0.,3^) 

PLOT (1.75,4. b*-3  ) 


, 14, 3H  X  , 14) 


R/S 

^-/S 

1  = 

J  = 

(I, 

(  I, 

Th 

XM 

XM 

1  = 

AM 

AM 

=  1 

J  = 

1  = 
CAL 
CAL 
.LT 

PL3 

32 

NUE 
TH 
1  = 

J  = 

CAL 
CAL 
.LT 
PLQ 

42 


CA 
CA 

I 

J 
J  ) 
J) 
E 

ES 
AX 

I 
AX 
IN 
./ 

J 

I 

E* 
E* 

•  "" 

1  ( 


L 
L 

0,  13 
1»  J3 

=  XO 
=  YO 
GRID 


+  . 5* ( A0(  1  )**2 
+A0( I) *( B0( J) 
CURVES  AROUND 


-(BC{ J  ) 
+  S0{ I  )  ) 
AIRFGIL 


♦  S0( I ) )**2  ) 


H (I  X 1 , 3  ) 


X  1 ,  I  X  2 

1  (XMESHd,  3),  XMAX) 

1 ( XMESH ( I, 3) >  XMIN  ) 

( xrAX-XMIN) 

1,  J3 

=  3 

0,  13 

(XMESh(  I, J  )-XMIN  ) 

(  rMESH(  I,  J  )-  Y:-lESh(  IG,  3)  ) 

3.75.DR.XP.GT.4.7b.CP .YP.LT.-4.5.DR.YP.GT.4.5)  GO  TO  33 

XP,  YP,KP  ) 

=  2 


E  GRID  CURVES  RADIATING  FRGM  AIRFOIL 

IU,I3 
=  3 

J1,J3 
E*  (XMESh( I, J)-XMIN )  ' 

E* (YMESH( I, J )-YMESH(I0,3) ) 
.-3.75.CJR.XP.Gr.4.75.QR.YP.LT.-4.5.0R.YP.GT.4,5) 

T(XP,YP,KP)  Mr... 

=  2 


=  3 


GO  TO  43 


NUE 

PLO 

SYM 

SYM 

PLQ 

FRA 

N 

PLOT (0.,0.,999) 


T(-1.75,-4.5',-3  ) 

B0L(-2.C,-1.50*.07,3,0.,-1) 

BuL(6.5,-1.50,.07,3,0.,-l) 

T(-2.5,-2.00,-3) 

ME(1) 
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RETURN 
EKO 


C 
C 
C 

c 


oUBRCUTINE  USTADI 
UNSTEADY  TRANSCMC  POT 
WITH  FIRST  QRDEP  kADlA 
SHEARED  PAKAdGLlC  CQOR 
DIRECTION  IMPLICIT  SCH 
GM(  122,  3b)  ,0, 
, AU( 132  ),A1( i 
33( 3b),H\,Hi 
AL>UTli^»CB*S 
DX,DY>uT>DXX 
RAD, PI>ALSiA 
, AMPLM»FRECRM>F ASAGh* 
,FREQR> IPSURE 


CGi^lMON/A/ 
1 

2  f 

3  > 
COMMGN/H/ 
CUMMON/J/ 

1 
2 


COMMON/K/  I0»  Il»  12*  I3» 

CUMMON/WAKE/  NIT>wG(13 

DIMENSION  C(132)»F(132 

COMPLEX   WA 

DDT=  l./DT 

DDXX=  l./DXX 

DDYY=  l./DYY 

NS         =  0 

AAO  =  1. 

RG  =  0. 

IG  =0 

JG         =0 

***** 

Y-SwEEP 
***♦♦ 

IM=  10  -1 

IMM=  10-2 

C  (  IM)=  0. 

C(  IMM)=  0. 

E  {  I M  )  =  0  . 

F(IMM)=0. 

F  (  I M  )  =  0 . 

FCIMM )=0. 

UPPER  BOUNDARY 

J  =  J3 

1=     10 

ANG  =  PI*.7t) 

Y»SO(I )+B0( J) 

X=     A0(I) 

Hh«     X*X    +     Y*Y 

DHH=     l./HH 

XT=    -.;?*Y*(ALT  +  CETAT  ) 

YT=    .5*X*(ALT+Ct TAT)    - 

H=SaRT(HH) 

Dh=    l./H 


ENTIAL    FLOW    EQUATION     IN    QUASILINEAR     FORM 
TIQN     aUUNDARY    CO^DITIONS     IN    MOVING 
DINATES     ARE     SOVLED    BY     AN     ALTERNATING 
EME     WITH    Y -SWEEP     FIRST 

(132*36),GN(13Z,36)*S0(132),S1(132)»S2(132) 
32)»A2(132)*A3(132)>bC(3b)»Bl(36)*i32(3b) 
,IX1,IX2>KSYM,FMACH, ALPHA, CA,SA,FM AC H2 
i3,NS,RC,IG,JG 

>dyy>dtt/dxy,0xt,oyt,tsr 

lt,altt,ampla»freqra,fasaga,fmachs,fmacht 

cetas,cetat,cetatt,amplc»freqrc,fasagc,ceta 

J1>J2,J3 

?) 

)  ,F(132) 


+     DriH*(CA*X    +     SA*Y) 
•DHH*  (CA*Y-SA*X  ) 
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GI 
GJ 
GX 
GY 
U  = 

V  = 
QQ 
CH 
Fl 
AA 
AA 
A  = 
WA 
U  = 

V  = 
AV 
TG 
TG 
IG 
TG 
Yl 

L 

dl 
DI 
tl 
CI 
AI 
GA 
bE 
AL 
C( 
El 
E( 
DO 
AN 

V  = 
X  = 
HH 
DH 
XT 
YT 
H  = 
DH 
GI 
GJ 
TG 
TG 
GX 
GY 
U  = 

V  = 
CO 
CH 
FI 


(G( I+l* J)-G( I, J ) )*2. 
(G( I, J)-G( I> J-1 ) )*2. 

=    Aid  )*    GI       -S1(I)*E1  (J  )*    GJ 
Bl( J )*    GJ 


GX*DH 
GY*DH 
=G*U+V* 
AIN=  XT 
T=  GM{  I 
«AA0-.2 
= AMAXl ( 
SQRT(AA 
=CMPLX( 
U+REAL 
V  +  AI 
=  \/-U*S 
1=  GI 
J=  GJ 
'M=    2  .  * 
MJ=  2.* 
=  -GM(  I 
-2 
=  0. 
=  0. 
'    0. 
=  DT*OH 
=  l.-CI 
MA=  DI 
OA=bI-C 
FA=  1./ 
I) =  (CI 
I)=  EI* 
!)=( YI- 
1  1=  I 
G=  .5*P 

so(  i)  +  b 

AO(I  ) 
=  X*X  + 
H=  l./H 
»  -.5*Y 
=  .5*X* 
SQRT(HH 
=  1  ./H 
=  &( I+l 
=  2.*(G 
J=  GJ 
MJ=  2.^' 


GX*DH 

GY*DH 

=U*U+V* 

AIN=  XT 

T=  GM(I 


V 

*GX  +  YT*GY 

»J)  ♦CUT  +  CHAIN 

♦  QQ-.'t  +  FIT 

AA> . 0001  ) 

) 

CQS{ANG)>SIN(ANG))*CMPLX(A,C.) 

(WA)   +H*XT 

MAG(WA)  +  H*YT 

1(1  ) 


(Gi'^l  I  +  l^  J)-Gi"l(  I,J  )  ) 

(Gf,(  I*  J)-Gh  (  I»  J-1)  ) 

»J)  +  DT*(U*TGMI*A1 ( I )+AV*TGMJ*Bl( J ) )*DH 

.*DT*(U*TGI*Al(I)+AV*TGJ*t]l(J))*DH 


*U*2.*A1(  I) 


( I-2)*GAMA 

(  AI-3EDA«-C(  I-l  )-GAMA*E(  1-2)  ) 

-BEUA+F  (  I-l  )  )*ALFA 

ALFA 

BEDA*F(I-1)-GAMA*F{I-2))*ALFA 

1>  12 

I 

(j(J  ) 


'■'  ;,  T  ;i 


H 

*(ALT+CETAT)  +  DHH*(CA*X  +  SA*Y) 

(ALT+CETAT)  -OHH*  (  C  A*  Y-S  A*X  )        • ': 

)  •;  t  ,  .  ♦ .  . 

>j)-G(l-i^j)  i;}  ^iL)?;! 

( I»J )-G( I^ J-1)  ) 

(GiM(  I,  J  )-&M(  I, J-1  )  ) 
=  Aid)*  GI   -S1(I)*B1(J)*  GJ     ^  ♦ 
=   Bl( J  )*  GJ 


V 

*GX  +  YT*GY 

*J)  *DDT  +  CHAIN 


(  '. 


I  '  o> 


H  1 1  »'  , 
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AA=4A0-. 2+QQ-.4+FIT 
AA=AMAX1 ( AA^.OOOl) 
A=S3kT( AA) 

WA  =  CMPLX(CQS( ANG)>SIN{ANG)  )*CMPLX(A,C.) 
U=    U+REAL(WA)       +H*XT 
V=     V    +     AIMAG(WA)     +    H*YT 
AV=     V-U*S1(I) 
IF(U.  LT.O.  )    GCJ    TC    2 
TGI=    2.*{G{ I> J )-G{I-l> J  )  ) 
TGMI=    2.*(GM(  I»J  )-GN.(  I-1>J  )  ) 
-    2.«A1(  I  ) *DT*OH*U 
-31 


BI  = 
AI=  1 
CI'  0 
GO  TIJ 
TGI  =  2 
TGMI  = 
CI  = 
AI  = 
BI  = 
YI  = 


*{G( I+l, J )-G( I, J ) ) 
2,*(GN(  I  +  i>  J  )-GM(  li  J ) ) 
2 .*A1  ( I )*DT*OH*U 

l.-CI 

0. 

-GM(I»J)     +    L)T*(U*TGMI*A1(  I  )+AV*TGf1J*31  (  J  )  )*DH 
1  -2.*DT*(U*TGI*Ai(l)+AV*TGJ*Bi(J))*UH 

DI=  0. 
tl=  G. 
GAMA=     DI 

bE0A=3I-C( 1-2 )*GAMA 

ALFA=     l./iAl-BEDA+ClI-D-GA-iA  +  Ed-d)) 
C(I)=     (CI-GtOA*E ( I-l ) )*ALFA 
E(I)=     EI*ALFA 

F(I)=(YI-BEDA*F(I-1)-GAMA*F(I-2))*ALFA 
.    CONTINUE 
1=     13 

ANG=     .25*PI 
Y=SO(  I)+BC( J) 
X=     AO  (  I  ) 
HH=     X*X     +    Y<=Y 
OHH=     l./HH 

XT=    -.5*Y*(ALT+CETAT)     +     DHH*(CA*X    +     SA*Y) 
YT=     .5*X*( ALT+CETAT)     -DHH* ( C A*Y -S A*X ) 
H=SQPT(HH) 
OH^     l./H 

GI=     (G(I»J )-G(I-l>J ) )*2. 
GJ=     (G(  I» J )-G{ If J-1 )  )*2. 

GX  =    Aid)*    CI       -S1(I)*B1(J)*    GJ 

GY  =       Bl( J )*     GJ 

U=     GX*DH 
V=    GY*DH 
QQ=     U*U+V*V 
CHAiN=     XK'GX     +    Yr*GY 
FIT=     GM( i> J )     ♦DDT    +    CHAIN 
AA=AA0-.2*QQ-.^*F IT 
AA=     AMAXK  AA>  .001  ) 
A=     5CRT(AA) 

WA  =  Ci^PLX(COS(ANG)^SIN(ANG))*CMPLX(A,0.) 
L)=    U+RFAL(WA)        +H*XT 
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V  = 
AV  = 
TGI 
TGJ 
TGM 
TGM 
YI  = 

CI  = 

Di  = 
tl  = 
bl  = 
AI  = 
GAM 
HED 
ALF 
C  (  I 
E  (  I 
F(  I 
CG  = 
CCG 
OJ 

r  . 

1  - 

DG  = 
CG  = 
CCG 
GN{ 
J  = 
LFF 
1  = 
ANG 
Y-b 
X  = 
HH  = 
DHH 
XT  = 
YT  = 
H=S 
DH  = 
GI  = 
GJ  = 
GX 
GY 
U  = 

V  = 
Q0  = 
ChA 
FIT 
AA  = 
AA  = 
A=S 
W  A  = 
U  = 

V  = 


V  f 

V-U 
=  CI 
=  GJ 
i=  2 
J=  2 

-GM 

0  . 
0. 

c . 
-OT* 

1.- 
A=  D 
A=3I 
A=  1 
)=  ( 
)=  E 
)  =  (Y 

U. 
=  0. 
H  K  = 
I3  +  I 

CG 

=  DG 

If  J) 

J-1 

T  BD 

TO 

=  P1 

D(  I  ) 

A0(  I 
X*'< 

=  1. 
-.5 
.5* 

QRT( 
1./ 
{G{ 

G(I» 


AIMAG( WA)  +  H+YT 
*Si  (  I) 


.*(GM(  1,J )-GM(  I-l^  J)  ) 

.»(GM(  I, J)-GM( I* J-1  )  ) 

(I/J)  ♦  CT*(U*To!"^I*Al(  I  )+AV*TGMJ*Bl(  J  )  )*DH 

-?.*DT*(U*TGI*A1(I)+AV*TGJ*31(J))*DH 


DH*U*2 .*A1 (  I  ) 

BI 

I 

-C  (  1-2 

./  {  AI-l 

C  I  -  3  E  L 

I*ALFA 

I-EEDA^FlI-D-GAMA^Fd-ZJJ  +  ALFA 


)  *GAMA 

PEDA*C ( I-i)-GAMA*E  (  1-2  )  ) 
:A*E(  I-l )  )*ALFA 


I0>13 
0-K 

F (I )-C (I )*CG-E ( I )*CCG 

=  CG 

U  N  D  A  R  Y 


+E0( J ) 
) 

+  Y*Y 
/Hh 

*Y*  {  AL 
X*( ALT 
HH) 
H 

I+l» J) 

J+1  )-& 

=  A 


GX*D 

GY*D 

U*U  + 

IN- 

=  6M 

AAO- 

AMAX 

QPT( 

CMPL 

U+RE 

V  + 


H 

H 

V*V 

XT*GX 

(I/J  ) 

.2*00- 

1  ( AA,  . 

AA  ) 

X(CLS( 

AL( WA) 

AlMAG( 


T+CETAT)  +  DHh*(CA*X  ♦  SA*Y) 
+CE1AT)  -DHH*(CA*Y-SA»X) 


-G( I, J  )  )*2. 

(  I>  J-1) 

i(  I)'*  GI   -Sl(  I  )*E1  (J  )*  GJ 

Bl( J)*  GJ 


+  YT*GY 

*DDT  ••■  CHAliN 

.^*FIT 

OOCl) 

ANG) >S IN( ANG)  ) *C MPL X ( A, 0 .  ) 

+  H*XT  -  ■   - 

WA)  +  H*YT 


'   '  ..4.   I 


■''  ■■  ,■;  '■ 


,  (^ 


!  i  »•  1^  V  ; 


,(i* 
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AV 
TG 
T& 
IF 
TG 
TG 
GO 
TG 
TG 
YI 
L 

BI 
DI 
EI 
CI 
Ai 
GA 
BE 
AL 
C  ( 
E( 
F( 
IN 
DO 
FX 
Y  = 
X  = 
HH 
DH 
XT 
YT 

cx 

CY 
XT 
XT 
YT 
YT 

cx 

CY 
BX 
BY 
XT 
L 
YT 

H  = 
DH 
GI 
GM 
GJ 
Gh 
GX 
GY 
U  = 


1=  GI 
MI=  2 
(  AV. 
J=  2. 
MJ=  2 
TD  7 
J  =  2.* 
MJ=  2 
=  -GM 

=  0. 
=  0. 
-  C. 
=  DT* 
=  1.- 
MA=  D 
DA=bI 
FA=  1 
I  )=  ( 
I)=  E 
I  )  =  (Y 
TFRIO 

8  1  = 
=  1.  +  S 
S0(  I) 

AO(I 
=  X*X 
H=  1. 
»  -.5 
=  .5* 
»  Y*Y 
=  2.* 
X=  (C 
Y=  -. 
X=  -X 
Y=  XT 
=  X** 
=  3.* 
=  2.* 
=  CA* 
T  =  -.5 
+  FMAC 
T=  .b 
♦  Fr'AC 
SQRT( 
=  1./ 
=  G(  I* 
1=  GM 
=  G(I» 
J=  GM 


=  V   -U*S1(  I) 

.*(GM( I+l» J)-GM( I, J) ) 
LT.O. )  CC  TO  6 
*(G( I* J)-G(I> J-1)  ) 
.*(&M(  I,J )-GM(I>J-l  )  ) 

{G( I>J+i)-G(I>J  )  ) 

.*(GM(  I, J  +  1)-GM(  If  J  )  ) 

(I*J)  +  DT*(U*TbMI*Al(I  )+AV*TGMJ*til(  J  )  )*DH 

-2.*CT*(U*TGI»A1(1)+AV*TGJ*B1(J)) *Dh 


OH* 

CI 

I 

-C  ( 

./  ( 

CI- 

I*A 

I-B 

R 

II 
1(  1 
+  B0 
) 

+ 
/HH 

*  Y* 
X*( 
-X* 
X*Y 
A*C 
5*  { 
TY 

X 

3-3 

X*X 

SA* 

CA- 

*  Y* 
HT* 
*X* 
HT* 
HH) 
H 

1,  J 
(1  + 
J+1 
(  I* 


U*2.*A1 ( i) 


1-2 )*GAfA 

Al-BEDA*C(i-i)-GAMA*E(I-2)) 

BEDA*F(  I-i)  )*ALFA 

LFA 

EDA*F(I-1)-GAMA*F(I-2))*ALFA 

,  12 

)**2 

(J) 

Y*Y 

(ALT+CETAT)  +  DHH*(CA*X  +  SA*Y) 
ALT+CETAT)  -uHH * ( C A ^Y-S A*X ) 
X 

X-CY*i)A)*DHH*DHH 

ALT+CETAT )-OHH*DHH* (SA*CX+C A* CY) 


.*X*Y*Y 

*Y-YV-*3 

CA 

SA*SA 

(ALTT+CETATT)-.2  2*X*(ALT+CETAT)**2 

OHH*(X*CB+Y*SB ) -AL T*DhH* { X +5 A-Y *C A ) -DHH**3* ( C X* B Y+ BX*C Y ) 

(ALTT+CETATT)-.2i>*Y*(ALT  +  CETAT)**2 

DHH* (X«SB-Y*Cd)+ALT*DHH*(X*C^+Y*SA)+DHH**3*(BY*CY-bX*CX) 


)-G(  I-1*J) 

ifj  )-GM(I-l, J) 

)-G( I,o-l) 

J  +  1  )-Gf''(I,  J-i) 

=  Aid  )*  GI   -SI  (1  )*[;1(J  )*  GJ 

=   Bl( J )*  GJ 


GX+DH 
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AU  = 
A\/ 
LR  = 
Vk  = 
AVR 
UUK 
VVR 
Q  u,  P 

s  = 
Ih( 

T=l 
IF  ( 

UU  = 

uv  = 
vv  = 

ChA 

FIT 

AA  = 

AA  = 

AB  = 

oil 

GIJ 

GJJ 

GMI 

GMI 

GMJ 

kOT 

CX  = 

CY  = 

AX  = 

AY  = 

A  k  = 

P  = 

IF( 

CfcN 

AXX 

AXY 

AYr 

YK  = 

AXY 

Yi^R 

BB- 

DI  = 

QI  = 

Al' 

CI  = 

EI  = 

GD 

TYP 

i^S 

K 

IM 


GY*DH 

L  +  V«S1  (  I  ) 

=  V   -b*Sl (I ) 

XT*H+  U 

YT*H  +  V 
=  VR-Uf^*Si  (  I) 
=  Uf<*UF 
=  VR*VR 
=  OUK  +  VVR 
1. 

Uk. LT.O. )  S=  -1 


AVR  .L 

iJ*U 

U*V 

llb'  +  VV 

IN=  X 

=  G,-^( 

AAO-. 

A  M  A  X  1 

Al(  I  ) 

=  (G(  I 

=  G(  1  + 

=  (G(  I 

I  =  (  G  I'l 

J=GM( 

J=  (GM 

ATEO 

XTT  + 

YIT 

Aid 

El(  J 

AX*G 

CQ  * 

OCR.G 

Tk  AL 

=  ( AA 

=  -2. 

=  B2( 

R  + 
=  -2. 
=  AXX* 
.5*0 
C. 

bB*( 
1.-2 
BB<-  { 
0. 
TC  10 
t  DEP 


T.O. )  T=  -1 . 


T*Gx 
W  J  ) 
2*G{J 
(  AA* 
*B1( 
+  1*  J 
1*  J  + 
*J+1 
(  I  +  l 
I  +  l» 
(I,  J 
CUQP 
2.* 
+  2. 
)*CX 
)«{C 
I  + 
(U*X 
E  .  AA 
DIFF 
-UUR 
*A6- 
J)  * 
AXX* 
*Uk* 
GMII 
TT*D 


+  YT*GY 

*l)DT  +  CHAIN 
-.^*FIT 
.0001 ) 

J ) 

)-2.*Gl  I*  J) +G(  I-li  J  )  )*[jDXX  +  A3(I)*GI 
l)-G(I  +  lfJ-l)-b(I-l>J  +  l)+G(I-l»vJ-l) 
)-2.*G( 1, J )+G( If J-i )  )*DUYY  +  b3(J)*GJ 

*  J  )-2.*Gh(  W  J) +Gri(I-l,  J)  )*DDXX  +  A3(I)*GMI 
J  +  l)-G^',(I+i*J-l)-G^(I-l*0  +  l)+GM(I-l,J-l) 

*  1)-2.*GM  (i,J  )-^Gr  (I  .  J-1)  )*DOYY  +  B3(J)*GMJ 
DILATES  TERM? 

(U*XTX+  V*XTY)*DH 

*  (U^YTX  +  V*YTY  )*[;H 


Y-CX*S1 {  1)  ) 

AY*Gj 

+V«Y)*DH  -( AA-OUR )*S2( I ) *GY  -HH*wR 

)  GD  TD  y 

ERENCING 

)  ♦A2{  I) 

(  AA*S1  (  I  )  +  UR'^^AVP  ) 

( AA*FX-AvR*AVR ) 

Gil  +  AXY*oIJ  +  AYY*GJJ 

AVR*Ab  ' 

+  AXY*GhIJ  +  AYY*GMJJ    '    '  ' 
HH*A2 ( I )♦ ( UUP-AA ) 


DDXX-A3( I ) ) 
.*B8*D0XX 
D0XX+A3( I ) ) 


ENDENT 
=  NS 

=  I   -K 


DIFFERENCNG 
+  1 
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12 


13 


IMi^  =     IM       -K 

L  =  T 

JM=J-L 
■JMrl  =  JM-L 

AUR=     UR+    VR*S1 ( I ) 
AQ=     AA/QQR 
bXX=    VVR*A2(I) 
BXY=     -2.*AB*V)<*AUP 
BYY=     AUR*AUR*B? ( J ) 
GNN=BXX*GII  +  BXY*olJ-»-3YY*GJJ 
GriNM=     BXX*GMI  I  +  BYY*G^  J  J 
GIJ^  =  G(  I*  J  )-G  (  I!^*  J)-G{  I  »  jM  )+G  (IM>  J,M  ) 
GMIJM=    GM(  U  J  l-GMC  IM>  J  J-G.'^  (  I^  J^  )  tGrK  IM,  JiM) 
IF (     JMM.GT  .  JJ  )     GG    T D    11 

bJJM=  (G(  I>  J  )-2.vG(  W  Jil) +G(  I,  JNf^  )  )*ODYY     +    B3(J)*GJ 
GMJJI«!={     GM  (  I»  J  )-2.*Gi1  I  I  »J.1) +GK(  I>  JMM  )  )*OOYY     *    B3(J)*GMJ 
GO    TD    12 
CJJM=GJJ 
GMJJM=     GMJJ 

IF{     IMM. LT.  IC  .Ok.  IHM.GT.  13 >     GO    TJ     13 
GIIM  =  {G(  I>  J  )-2.*G(  IM*  J  )+G(  IMi'-l,  J  )  )*DDXX     +    A3(I)*GI 
GMIIM=     (Gh(I>  J)-2.=i--GM{  IM^J  )+GM(  IMM,  J  )  )«OOXX     +     A3{I)*GMI 
GO    TD    14 
GIIM=     Gil 
GMIIM=    GMll 
AXX=     UUR*A2( I ) 
AXY=     6.  *S*T*lJR*AVR*Aii 
AYY=     AVR*AVR*B2(J) 
GSS=AXX*GlIM+AXY^blJN  +  AYY=i 
GMSS=     AXX*GM,II,'*i+AYY*&KJ  JM 
Yk  =     ( AQ       -1. )*jSS 

YMR  =  Aa*  (G^S.S  +  G,^,iNlN  ) 

GMSS  =     AXX*Gl'ilII'.  +  AXr♦GMIJ^'  +  AYY*GMJJM 

YMR=     YMR    -Gi^SS 

3d=     .5*DTT*0hH*UUF*(l.-Au)*A?(I) 
CC=     -.5*D1T*0HH*AG*VVR*A2(  I) 
BOCC=    BB+CC 
1F(     UR.LT.O.)     GLi    TU     13 

l.cQ.  ID     GU    TO    Ifa 

B  B  *  D  I)  X  X 

D0XX*(CC-?.*[3b)     -A3(I)*bbCC 

1.    +DCXX* ( 3B-?.*CC  ) 

CC*DDXX    +    A3( 1 )*BBCC 

0. 
TO    10 


ib 


G  J  o 


+AG*GNH        +R 


16 


IF  ( 

DI  = 
Bl  = 
AI* 
CI  = 
EI  = 
GO 
IF( 
DI  = 
BI^ 
AI  = 
CI  = 
tl  = 
GU 
31  = 
31  = 


i.co.12)  g;:  to  i6 

c. 

CC+UDXX-A3 ( I ) ♦BBCC 
l.+DDXX*(B8-2.*CC ) 
DDXX*(CC-2.*36)  +  A3(I)*bBCC 
BR*DDXX 
TO  10 
0. 
B3CC*(0DXX-A3(  I  )  ) 
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ai=     1.     -2.*DDXX*BBCC 

CI=     bBCC*(nDXX+A3(I)) 

EI=    0. 
C  ADVECTION     TEkMS 

C  UPWIND    DIFFFr<tNCIN& 

10    Yr=AVR*a.*31(  J)*(G^1(  I,  J  )-GF:(IiJ-l)  ) 

IF(AVP.LT.O.)     YY=AV3*31(J)*(GM(I>J+1)-GM(I*J))*2. 

iJti  =  Uk*DT*Lh*2  .'■^Al  (  I  ) 

IKUP.LT.O.)     GU    TC    18 

yy=     YY+UK*2.*A1(I)*(GM(I,J)-GM(I-1,J)) 

'31=     bl-BB 

AI=     AI+S3 

GLJ    TG    19 

18  YY=     YY  +  UR*2.*  All  1)»  (CI1{  I  +  i>  J  )-bM  (  1,  J  )  ) 
CI=    C1+3B 

AI=     AI-3B 

19  ri=     GM(  I*  J  )+OTT*(  Y'R-.  5*Yf'R  )  ♦DHH    -DT*YY*DH 
GAl^A=     DI 

EEDA=BI-C (  I-?)*GAMA 

ALFA=     l./(Al-bEUA*C(I-i)-GA,^A*E(l-2)) 
C{I)=     (CI-tiEDA*=L  (  I-l  )  )*ALFA 
Ed  )=     EI*ALFA 

F(l)  =  (YI-BEDA*F(I-l)-GAi>lA*F(I-2))*ALFA 
8     CONTINUE 
C  RIGHT    Bau^sDAkY 

1=     13 
A  N  G  =     0  . 
Y=SO(  I) +B0( J) 
X=     AC)  (I) 
HH«     X*X     +    Y*Y 
DliH=     l./HH 

XT=  -.5)*Y*  (ALT  +  CFTAT)  +  UHH*(CA*X  +  SA*Y) 
YT=  .^*X* ( ALT+CET AT)  -DHK* { C A* Y-SA*X ) 
h=SQFT( HH) 
Ch=  l./H 

GI=  (G(I* J  )-G(  I-1,J )  )*2. 
GJ  =  G(  I^  J  +  1  )-G  (  I* J-1  ) 

GX         =  Aid)*  GI   -Sld)*El(J)*  GJ 
GY         =   31(J)*  GJ 
V=     GX*DH 
V=  GY*OH 
GQ=  U*U+V*V 
CHA1N=  XT*0X  +  YT+GY 
FI T=  GM{  I> J  )  *DDT  +  CHAIN 
AA=AA0-.2*QQ-.^*FIT 
AA=  AfiAXl  (  AA..0001  ) 
A=SQRT( AA) 

■wA  =  CMPLX(COS(  ANG  )*S1N(ANG))*CMPLX(A,0.) 
U=  U+REAL(WA)   +H*XT 
V=  V  +  AIMAGCWA)  +  H*YT 
AV  =  V   -U*S1(I) 

TGI=  GI 

TGMI=  2.*(GM{  I*  J  )-GMd-l,  J)  )     ,  ,   ,       (•!.  .  ] 
IF (  AV.LT.O. )  GO  TO  2  0 
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20 


TGJ=  2.*(G(  I, J)-C(I» J-1 ) ) 
TGMJ=  2.* ( GM( I> J)-GM( I^ J-1) ) 
GG  TC  21 

TGJ=2.*(G(I.J+1)-G(I>J)) 
TG.^J=  2.*(Gh(  I»  J+1)-GM(I»J  )  ) 


CG 
3  ) 


GO  TO  b 


21  YI=  -Gld^J)  +  DT*(U*TGMI*  Al  (  I  ) +AV*TGMJ*dl  (  J  )  )*0H 
1         -2.*UT*(U*TGI*A1(I)+AV*TGJ*B1(J))*DH 

DI=  C. 

cl=  0. 

CI=  0. 

BI=  -0T*DH*U*2.*A1(  I  ) 

AI=  l.-BI 

GA'^A=  DI 

BtDA=BI-C ( 1-2 )*GAMA 

ALFA=  1. / (AI-SEDA*C ( I-l )-GAMA^e  {  1-2)  ) 

C(I)=  (C1-BEDA*E( 1-1) )*ALFA 

FU)=  EI*ALFA 

F(I)=(YI-bE0A*F(I-l)-GAMA*F(I-2))*ALFA 

Co=  0. 

CGG=  0. 

DG  22  K  =  10>  13 

1=  I3+I0-K 

DG=  CG 

CG=       F  (  I)-C  (  D-i-CG-t;  (  I  )*CCG 

CCG=  DG 

22  GM(i>J)= 
IF(  J.GT 

C      ***** 

C      X-SwbEP 
c  *  *  *  *  * 

C(l)  = 
C(2)  = 
E  ( 1  )  = 
E(2)  = 
F{1)  = 
F  (2)  = 
C       LEFT  BOUNDARY 
1=  10 
DG  23  J= 
I  F  (  J  .  F  Q 
iF(  J  .EQ 
AK6=  PI 

3J=  G(I» J  +  l)-G(  I* J-1 ) 
GO  TC  26 
ANG=  1.2  5*PI 
GJ  =  G(I*  J  +  1  )-G(  1>  J-1  ) 
GO  TO  2b 
ANG=  ,7&*PI 
GJ=  2.*(o( I»J  )-G(  I^ J- 
Y  =  SO(  I  )  +  B0 (J) 
X=  AO(I) 
HH=  X*X  +  Y*Y 
UHH=  l./HH 
XT=  -.5*Y* ( AL T+CFTAT) 


0, 

0, 

0, 

0 

0, 

0. 


Jl»  J3 
Jl)  GO 
J3)  GO 


TO 

TO 


25 


2^ 


25 


?6 


+  UHh*(CA*X  ♦  SA*Y) 
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25 


27 


2Q 


Z3 


30 


YT 

Mli 
H  = 
DH  = 
GI  = 
GX 
GY 
U  = 

v= 

00  = 
ChA 
FIT 
AA  = 

AA  = 
A=S 

f,^  = 
u  = 
v  = 

AV  = 
U  ( 

IF{ 

PI  = 
AI  = 
Cl  = 
GO 
CI  = 
AI  = 
BI  = 
Yl  = 
Dl  = 
EI  = 


^A0( 
,QPT 
=  1  . 
■■     (G 


I  ) 

(H 
/H 
(I 


♦(ALF+CETAT)  -OHH* ( C A* Y -SA* X  ) 

*A0(  I  )  +  Y* Y 

H) 

+1»J )-G( I> J ) )*2. 

=  Al( I ) *  Gl   -SI ( I ) *B1 ( J )*  GJ 
=   Bl( J )*  GJ 


BE 
AL 
C( 
E( 
F  ( 
CO 

cc 

CG 
DfJ 
J  = 
IjG 
CG 
CC 
GN 
IN 
uO 
FX 

eu 

IF 

Y  = 
X  = 


GX* 
GY* 
U*U 
IN  = 
=  G 
A  AG 
AM  A 
UFT 

c^'iP 
u+s 

V    + 

V- 

J. 

J. 

AV 

1. 

G. 
TC 
DT 
1  . 

c. 

GN 
U. 
0. 

A  = 

A  =  3 

A=: 
)  = 
)  = 
)=    ( 

TIN 
=    0 

0. 
3  0 
J3 

CG 

F( 
=  D 
If  J 
ERI 
31 

1. 
32 

J  . 
0(  I 
A0( 


DH 

DH 

+  \/ 

X 

XI 
(A 
LX 
EA 
A 
U* 
£Q 
EQ 
.L 
OT 


*V 

T*G 

If  J 

2*0 

(  AA 

A) 

(C3 

L  {A 

IMA 

Si  ( 

.Jl 

.  J3 

T.O 

*DH 

I 


X     t     YT*GY 

)     *LiUl     +     CHAIN 

0-. 4*FIT 

,. pGOl) 

S(ANG)fSIN(ANG))*CriFLX(AfO.) 

A)        ■H-*XT 

G(  *<A)     *■    H*YT 

I  ) 

)     GU    TO    27 

)    Gc  TG  ze 

.  )     GU    TD    27 
*A\/*2.*tJl  (J  ) 


29 

*0H*AV*2.*B1( J ) 
-CI 

(  I»  J  ) 


DI 

I-C ( J-2  )*GAMA 

l./(AI-eEDA*C(J-l)-GAMA*E(J-2)) 

(C  I-i3E0  A*E  (  J-1  )  )  *ALFA 

E  I  *  A  L  F  A 

YI-BE0A*F(J-l)-GAf^A*F(J-2))*ALFA 

UE 


K=     Jlf J3 

+     Jl     -K 


J  ) 
G 

)  = 
OK 
1  = 
+ 
J  = 
EQ 
)  + 
I  ) 


-C(  J)*CG-E(J  )*:CG 

CG 

Ii>  12 
SKI)  **2 
Jl>  J3 

.J3)    GJ    ra   33 

B0(  J  ) 


■  I 


♦  < 
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3^ 


-SI  (I )*ei(J )*    GJ 


hH=     X*X     +     Y*Y 

OHH=     l./HH 

XT=    -.&*Y* (ALT  +  CETAT  )     +     DHH*(CA»X    +     SA*Y) 

YT=     .5*X*(ALT+CETAT  )     -DhH* ( C A* Y-S A* X  ) 

H=SQRT(HH) 

0H=    l./H 

GI=G(I+1*J )-G( I-l^ J ) 

GJ  =  G(  I» J+1  )-G(  I» J-1 ) 

GX  =    Aid  )*    GI 

GY  =       Bl(J)*     GJ 

U=     GX*DH 

V=    GY*DH 

AU=    U+V*S1(I) 

AV  =    V       -L'*S1(  1) 

URx     XT*H+    U 

VR=     YT^H    +    V 

AVk=    VR-UR*S1(I) 

UUR=    UR*UR 

\/VR=    VP*VR 

Q(jR=    UUR    +    Vyp 

CHAIN=     XT*GX    +    YT*GY 

FIT=     GX(I^J)     ♦our    +     CHAIN 

AA=AA0-.2*UQ-. ^*F I T 

AA  =  Ai^AXH  AA^.OGCi) 

IF(QQR.GE.AA)     GC    TG    3^ 

BB=     .b*DTT*DHH«(AVR*AVP-AA*FX)*B2(J) 

DI«    0. 

61=    bB*( DDYY-a3( J ) ) 

AI=     i.-2.*BB*DDYY 

Cl=    PB*  (DDYY    +     t>3  (  J  )  ) 

EI=>    C. 

GC    TO    36 

A(^=     ;iA/UQR 

AUR=    UR+VR*S1(I) 

BB=     .5*DTT*UHH*AVR*AVk*(1.-AQ)*62(J) 

CC=     -.5*DTT*DHM*AC*AUR*AUR*B2(J) 

t)BCC=    BB+    CC 

IF  (     J.Ea.'t  )     GG    Tu     38 

IF  (     J.FQ.J2)     GC    TO    39 

IF(AVR.LT.0.  )     GL.    TO    'tU 

0I«    BB+DOYY 

DDYY*  (CC-2.+dii)     -Li3(  J  )*B[CC 

1.     +    DDYr*(BB-2.*CC ) 

CC*DOYY    +     B3( J  )  ♦cBCC 

C. 
TC    36 


^0 


0l  = 

AI  = 

CI  = 

EI  = 

GO 

DI  = 

ril  = 

AI  = 

CI  = 

EI  = 

Gb 

0. 

CDYY*CC    -B3( J )*BBCC 
1.     +    DaYY*(BB-2.*CC) 
CDYY*(CC-2.*BE)+B3(J)*3BCC 
BB*ODYY 
TO    35 


38    DI»     0. 
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^1 


39 


^2 


36 


A3 


33 


IFi 
EI  = 
Bl  = 
AI^ 
CI  = 
&C 
Bl  = 

AI: 

CI  = 
t  1  = 

GC 
II-- 

IF 
01  = 
61  = 
AI  = 
CI  = 
Gu 
01  = 
31  = 
AI  = 
C  i  = 
A0( 
3b  = 
IF 
BI  = 
AI  = 
Gu 
CI  = 
AI  = 
GQ 
Ah; 

Y  =  : 
x  = 

HH  = 
DHf 
XT  = 
YT  = 

H  =  : 

DH  = 
GI  = 
GJ  = 
GX 
GY 
U  = 

v= 

u  ^- 
CHi 
FT 
AA  = 
AA  = 

A=: 

WA  = 
U  = 


AVR 

■  C. 

■  &B 
=  1. 

:   B8 

TG 

■  CC 
^  1. 

B6 

TG 
^  0. 

AVR 
:  F.P 
^  DD 

■  1. 
^  DJ 

TC 

■  G. 

^  53 
^  1. 

■  Bb 
'ECT 
^  DT 

AVR 
:  BI 
:  AI 

TO 
■■    CI 

■  AI 
TC 

p  -     • 

,0  {  I 

AG{ 

X* 

\'    1 

"■  • 

^  .5 
iJRT 

1. 
^  G( 

Z  • 


.LT.O.)  GG  TQ  '♦I 

CC*(D0YY-B3 ( J ) ) 

-2  .*E6CC*0DYY 
CC=f(  CDYY  ♦  B3(  J  )  ) 
35 
*DCYY  -33(J  )*BBCC 

+  LOYY*(  Bci-2.*CC  ) 
YY*  (CC-?..*BB)  +  •J3(J)*6bCC 
*ODYY 
36 

.LT.O. )  GU  TG  A 2 

*DDYY 

YY*(CC-2.*E3)  -33(J)*o&CC 

+ODV Y*(Ba-2.*CC) 

YY=i'CO  +  [^3(J)=f=BD:C 

36 


CC*  ( 
-2.* 
CC*( 

ICN 
*DH* 
.LT  . 
-BP 
+  B6 
46 

+  Bb 
-6'6 
46 

5*Pi 
)  +  B0 
I  ) 
X  + 
.  /HH 
b*Y* 
=^X*( 

(riH) 

/H 
I  +  l^ 
*  (G( 


DDirY-B3(  J  )  ) 
PBCC*Ll)YY 
DDrY  +  B3(  J)  ) 
TERMS 

AVR*2.*B1( J ) 
0. )  Gu  TG  43 


GX* 
GY* 
U*U 
IN  = 
=  G 
AAO 
AMA 
QPT 
CMP 
U  +  R 


DH 
DH 

+  V*V 
XT* 
M(I, 
-.2* 
XI  (A 
(  AA) 
LX(C 
EAL( 


(J  ) 

Y*Y 

(ALT+CETAT)  +  DHH*(CA*X  +  SA*Y) 
ALT+CtTAT)  -OHH*{CA*Y-SA*X ) 


J  )  -  G  (  I  - 1  ^  J  ) 

I^  J  )-G  (  I,  J-1)  ) 

=  Al(  I  )*  GI   -:ii(I  )*B1(J  )*  GJ 

=   B 1  (  J  )  *  G  J 

,  A  ^'  A  ..  -  Y  -'  ':   .  '  ■■>-;;:■ 

GX  +  YT*GY 

J  )  *DDT  +  CHAIN  '^  •  : 

C  C-  .  4  »  F  I  T  !  ■  I  •"  .  J  ■  ' 

A,.  OCOl )  ■    ' 

QS( ANG) »S IN{ ANG) ) *C MPL X ( A, 0 . ) 
WA)   +h*XT 


\  ■"  .  ' 


'  K 
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+  H*YT 


^8 


Af) 


32 


^9 

31 


51 


bZ 


'j2, 


V  +  AIMAG(WA) 
V-U*S1(I ) 

-DT*DH*AV*B1 ( J) fZ. 
l.-BI 


1  = 

1  = 


YI-BI*&M( 1,2) 
YI-DI*GM(I,2 ) 


-DI*GM(I,1) 


C. 

0. 

0. 

CN(I> J  ) 

J.EQ.3) 

J.EO. ^) 

DI 

I-C ( J-2)*GAMA 

1./  (AI-BEDA*C ( J-l>-GAMA*t  ( J-2 )  ) 

(CI-BfcDA*F ( J-1 )  )*ALFA 

E  I*ALFA 

)-GAMA*F (J-2 )  )*ALFA 


)*:CG 


V 

AV  = 

BI  = 

AI  = 

CI  = 

Dl  = 

EI  = 

YI  = 

IF( 

IF{  J 

GAMA  = 

BEDA  = 

ALFA  = 

C(J  )  = 

E(  J)  = 

F(J )=(YI-BEDA*r I J-) 

CONTINUE 

CCG=  0. 

CG=  0. 

DD  ^9  K=  Jl, J3 

J=  J3+J1-K 

D6=  CG 

CG=  F{J  )-C ( J) *CG-t ( 

CCG=  DG 

G  N  (  I  >  J  )  =  C  G 

CONTINUE 

RIGHT  BOUNDARY 

1=  13 

DO  bO    J=  Ji>J3 

IF(  J.EQ.Jl)  GO 

IF  {  J  .E3.J3)  GO 

ANG=  0. 

GJ=  G(I,J+1)-G(I, J-1) 

GO  TO  53 

ANG=  -.25*PI 

6J=  G( I» J+l)-G( I> J-1) 

GO  TO  53 

ANG=  .25*PI 

GJ=  2.*(G(  I, J  )-C ( I, J-1)  ) 

Y  =  SO(  I)  +  60( J) 

Y=SO(I )+B0( J ) 

X=  A0(  I) 

HH=  X*X  +  Y*Y 

UhH=  l./Hhi 

XT=  -. 5*Y*(ALT+CETAT)  +  DHH*(CA*X  +  SA*Y) 

YT=  .5*X* ( ALT+CETAT )  - JHH* ( C A* Y- SA*X ) 

H=SQRT(HH) 

DH=  l./H 

GI=  2.*(G(I,  J  )-G(  I-l, J  )  ) 

GX         =  Al(I)*  GI   -SKI  )*B1(J  )*  GJ 

GY         =   Bl{ J)*  GJ 

U=  GX*OH 

V=  GY*DH 

QO  =  U*U  +  V*\/ 


rc 

TO 


51 
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54 


56 


bO 


57 


58 


CHAIN 
FiT  = 
AA=AA 

AA=AM 
A=SQK 
w  A  =  C  M 
u=  U  + 
V=     V 


AV  = 
IH  ( 
1F( 
IF( 

BI  = 
Ai  = 
Cl  = 


V 
J 
J 
A 

1 

G 

TO 


=     xr*GX     +    Yr*GY 

Gi^d^J)     *DL/T     +     CHAIN 

0-.2*QG- . ^♦F IT 

AXKAA^.OOCl) 

T{  AA) 

PLX (CDS(  ANC) ^3IN( AN&)  )«CMPLX( A*0.  ) 

SEAL(WA)        +H*XT 

+     AIMAG( WA )     +     H*YT 

-U*S1  (I) 

.EQ  .  Jl  )     GO     TO     5A 

.EG.  J  3)     Gl?    TG    5  5 

V.LT.O.  )     GU    TO    5A 

-0T*DH*AV*2:.*B1  (  J  ) 

.-BI 


50 


CI  = 
AI  = 

ai  = 

YI  = 

Lil  = 

EI  = 
GAKA  = 
ri  t  0  A  = 
ALFA  = 
C(  J)  = 
E(  J  )  = 
F  (  J  )  = 
CuNTI 
CCG  = 
CG=  0 
DU  57 
J=  J3 
CG=  C 
CG=  F 
CCG  = 
GN(  I, 
UPOAT 
UC  58 
DD  56 
CG=  G 
G(  I>J 
GM(  I^ 
IF(  A 
RG=  A 
IG=  I 
JG=  J 
CGNTl 
I  F  (  N  I 
UTIM  = 
F  A  S  A  G 
ALPHA 
AL=  A 
ALT  = 


T*DH*AV*?.*H1 (J ) 
.-CI 


D 

1 

0. 

G  N  (  I  ^  J  ) 

0. 

0. 


DI 
31- 

1  . 

(C 

t  I 
(YI 
NUF 


K  = 

+ 

G 

{J  ) 

DG 

J)  = 

E  N 
1  = 
J" 

N(  I 

)  = 

J)  = 

bS( 

as( 


C ( J-?  )*GAhA 

/(AI-RECA*C(J-l)-GAMA*E(J-2)) 
I-QEDA*t ( J-i ) ) *ALF A 
♦  ALFA 
-dEDA*F(J-1)-&Ai'1A*F(J-?))*ALFA 


J1*J3 

Jl  -K 

-C( J) *CG-F( J  )*CCG 

CG 
EXI  RUN  DATA 

10,  13 

J1,J3 
>J  ) 
G(  I,  J  )+CG  ....-■',   -t  ^;  V 

CG 
CG).LE  .KG)  GQ  TO  58  .  i'  (.:  I 
CG) 


)  ■ 


NUE 

T.NE.O)  GO  TG  59 

UTIM  +  DT 

A=  FkEQRA*UTI.^ 

=  ALS  +  Ar^FLA*S  IN(  FAS  AGA) /RAD 

LPHA*RAD 

A'lPLA^FREwkA^COljiFASAGAj/kAO 


^  1  ,,  . 
\  f  '  .  ..  I 


196 


ALTT=  -AMPLA*Fi^EQRA*-«2*SI!>l(FASAGA)/RAD 

FASAGM=  FPEQRrt*UTir 

FMACH=     FMACHS     +    A  f.P  L  h*i  I  N  (  F  AS  AGM  ) 

FMAChT=     AMPLM*FI<buRM*CGS(FASAGr) 

FASAGC=  FkEQPC^UTIM 

CETARD=CETAS  +  A,";  PL  C  *S  i  N  (  F  AS  A  GC  )  /  R  AD 

CETA=  CETARO  ♦RA') 

CETAT=  AMPLC*FPEaRC*CQS(FASAGC)/RAL 

CETATT=  -AMPLC*FREGRC**2*5  1N(FASAGC)/RAD 

Cb=  COS(ALPHA) 

Sd=  SliM(ALPHA) 

CA=  FMACH*CB 

SA=  FMACH+SB 

FMACH2=  FMACH^+'Z 

WA<E  CONDITION 


i) 


■,L    TG  CO 


13 

3)-G(  IXI, 


59  IF(  NIT.EQ 
Gu  TO  61 

60  DO  62  1=  1X2^ 

62  WG{I)=  GdXZ, 

61  1=  1X2 
^G(I)=  WG(I) 

63  1=  I+l 

Y=SO( I)+B0(3) 

X=  A0(  I  ) 

HH=  X*X  +  Yi^Y 

DHH=  l./HH 

XT=  -.5*Y* ( ALT  +  Ce TaT  ) 

H=SQRT (HH) 


3) 


+  GiM(  1X2^  3)-GN(  1X1  >3  ) 


+  OHH*{CA*X  i-  SA*Y) 


YP  = 

HP  = 
GI  = 
IF{ 
GJ  = 

UP 


Y 

h 


3)  ) 


-SI  (I  )*bi(3)<'  GJ  )  /H 


G(I+1>3)-G( I-l»3) 

1.E0.I3)  GI=  2. *(G( If 3)-G(I-lf 

2.* {G( I*^ )-G( I»3)  ) 
=(Al( i)*  GI 
M=  NX<-  ^  -I 
Y=S0(M)+B0(3) 
Hh=AO(M  )*AO(:'i  )+Y-Y 
t<=SQRT(  MH) 
HN=  H 
Yi"l«  Y 

GI=  G(M+1, 3)-G(M-l>3) 
IF (  l.EQ.  13)  GI=  2.*( 
GJ=2.*  (G(h,4)-b(yi,  3)  ) 
U  .-^  =  (  A  1  (  M  )  «  G  I 

Y=  .5*(YP-Yrt) 
U=  .5*(UP-UM) 
h=  .5*(HP+HM) 
&F=  2.*DT*A1(  I  )*(U/H  +XT) 
WG { I  )  =  (  W  G (  I  ) +  B  F  *  WG ( I -i )  ) /  ( 1 .  +  B  f ) 
1F(  I.LT.  13)  GO  TO  63 
DLl  o7  1=  10,13 
CCG=  G(I,1  ) 


(M<-I,3)-C(I1,3)  ) 
-SKM)  »C'l  (3)*  GJ)/H 


CG- 

M  = 


G(I,2) 

N  X  +  ^  -  1 


197 


bb 


Ob 


b'f 


b7 


Ih  ( I .GT. IXZ)  Gu  in  bb 

IF  {  I  .LT.IXl  )  GC:  T[l  bb 

TANGtNTIAL  BGUNUARY  CUNDITIUN 

Y=SO( I )+BC(3) 

X  =  A  0  (  I  ) 

HH=  x*x  +  Y*y 

DHH=  l./HH 

XT=  -.5*Y*( ALT  +  C£TAT  )  +  DHH*(CA*X  +  SA*Y) 

YT=  .5*X*( ALT+CfcTAT  )  -OHH * ( C A *Y-SA*X  ) 

VBN=HH* ( XT*Si (  I  )  -YT) 

GI=  &(I+li3)-G( I-l^ 3) 

FX  =  1.  +  S1 ( i  )**2 

Bl S=  FX*ei ( 3) 

GXSXVB=  41 ( 1 ) *&I*S1( I )+VBN 

G(I»Z)  =  &(I>^)  -GXSXVa/BlS 

G(I*i)=  G(I*b)  -2.*GXbXVd/61S 

GO  TO  6^ 

G(I,2)=  G(M,^) +Wo( I  ) 

G(  I»  1 )=  G(M,5  )+Wol  I  ) 

GO  TO  64 

G(I*<^)=  GiM*'^  )-UG(M) 

G(I>1)=  G('^*5-)-wG(M) 

bM(  I, 2)=  G(  I»2)-CG 

GM(I»1)=  G(I>1)-CCG 

CuNTlNUE 

R  ETURN 

END 
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